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FOREWORD 

The  work  described  in  this  report  was  performed  in  the  Space  and  Surface  Systems  Division 
with  partial  support  from  the  Computer  and  Information  Systems  Division  of  the  Strategic  Systems 
Department.  Its  purpose  was  to  develop  a  new  algorithm  for  the  elIip>soidal  coverage  function,  and  to 
design  associated  Fortran  software  which  is  suitable  for  inclusion  in  a  high  quality  mathematics  and/or 
statistics  subroutine  library. 

This  document  was  administratively  reviewed  by  J.  L.  Sloop,  Head,  Space  and  Surface  Systems 
Division.  The  flowchart  on  page  13  was  prepared  by  Dottie  J.  Burgess  on  a  Macintosh  IIx  personal 
computer. 
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I.  INTRODUCTION 


Three  dimensional  spherical  coverage  problems  often  appear  in  the  study  of  weapon  evaluations 
for  aerial  and  submarine  targets.  Such  studies  require  the  capability  to  compute  the  ellipsoidal  coverage 
function  and  its  inverse.  The  ellipsoidal  coverage  function  P(R,  H,K,L,  u,  v,  w)  represents  the 
probability  of  an  event  occurring  in  Oxjy^z^ -space  within  a  sphere  SP  with  radius  R  and  center 
(H,  K,  L),  where 

SP:  (xi-H)2-»-(y,-K)2-j-(z,-L)2  =  R2  . 

The  random  event  occurs  under  an  uncorrelated  trivariate  normal  distribution  with  mean  (0,0,0)  and 
standard  deviations  u,  v,  w  in  the  Xj,  yj,  and  Zj  directions,  resjjectively.  Thus 


P(R,H,K,L,u,v,w)=  I  1 1  ir(x„  yj,Zj,u,v,w)  dxj  dyj  dzj. 


(1) 


where 


and 


Volume  of  SP 


3'(Xi,yj,Zj  ,u,v,w)  = 


(•^)^  u  V  w 


E[xi/(>[2  u)]  E[y,/(>f2  v)]  E[z,/(^^2  w)]  ,  E(z)  =  e  "  ^  (2) 


P(R,  H,K,L,u,v,w)  =  P(R,  K,H,L,v,u,w)  =  P(R,  H,L,K,u,w,v)  =  P(R,|H|,|K|,|L|,u,v,w).  (3) 
Setting  Xj  =  •'(2  u  X,  yj  =  -^  v  y,  Zj  =  ^2  w  z,  one  obtains 


L-l-R  K-l-Y  H-l-X 

P  =  E(z)  dz  I  E(y)  dy  |  E(x)  dx  ,  (4) 

■^L-R  "^K-Y  "^H-X 

where 

X  =^R  2_(L-V2wz)2-(K->f2vy)Y(^u)  ,  Y  =  J  R  *  -  (L  -  ^2  wz)  ^(>(2  v)  , 

and 

J  H  =  H/(a(2u),  K  =  K/(>|2v),  L  =  L/(^w) 

I  Ri  =  R/(>f2u),  R2  =  R/(V2v),  R3  =  R/(>j2  w)=:R. 

It  is  easy  to  show,  by  using  (4)  with  normalizations,  that  P  is  a  function  of  six  independent  variables. 

Rather  than  carry  out  the  numerical  triple  integration  of  (4)  directly  we  make  use  of  an 
available  computer  program  PKILL  (or  CIRCV)  that  yields  the  probability  Pc(r,  H,  K,  u,v)  of  an 
event  occurring  under  a  bivariate  normal  distribution  inside  a  circle  in  the  Oxjyj  -  plane,  with  center 
(H,  K)  and  radius  r.  P,.  is  the  two-dimensional  analog  of  P  [3,4,5, 6, 7]. 

Geometrically,  one  observes  that  P  can  be  obtained  by  considering  circular  slices  of  SP  parallel 
to  the  Oxy-plane.  For  a  fixed  z  in  [L  —  R,  L  -1-  R),  the  xy-integration  over  a  slice  of  radius  r  yields  xP^. 
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Weighting  with  E(z)/-^  and  integrating  the  result  over  z  in  [L  —  R,  L  -|-  R],  gives  P,  i.e., 

L4-R 

P=:i[  E(z)  P<.(r,H,K,u,v)dz,  r  =  (g) 

Another  useful  form  for  P  can  be  obtained  by  splitting  the  integral  in  (6)  into  two  pieces.  One  carries 
the  integration  from  L  —  R  to  L  and  the  other  from  L  to  L  +  R.  For  the  second  let  z  =  2L  —  t,  then 
combining  the  results  gives 

_ 

P  =  i  [E(t)-l-E(2L-t)]  Pc(r,H,K,u,v)dt,  r  =  ^R^ - (L -aJz w t) ^  .  (7) 

L-R 

The  symmetry  prop>erties  of  P  indicated  by  (1)  and  (3)  allow  H,  K,  L  to  be  taken  nonnegative. 
Since  the  integrand  of  P  is  positive  and  bounded,  the  order  of  integration  in  (1)  is  immaterial.  Thus,  as 
long  as  H  is  associated  with  u,  K  with  v,  and  L  with  w,  it  does  not  matter  which  is  called,  say,  L  and 
w.  For  example,  if  the  order  of  integration  is  chosen  so  that  the  original  x-integration  is  performed  last, 
then  if  initially  H  =  10,  u  =  5,  L  =  20,  w  =  7,  we  simply  let  H  =  20,  u  =  7,  L  =  10,  =  5.  In  this  way, 

we  can  refer  to  (6)  or  (7)  as  the  basic  representations  for  P,  where  L  and  w  are  always  associated  with 
the  z-integration,  with  the  understanding  that  the  original  order  of  integration  may  have  been  changed 
and  the  variables  H,  K,  L,  u,  v,  w  renamed  as  above. 

The  objective  in  this  report  is  to  expand  on  the  work  described  in  [8].  In  [8]  for 
-t-  4-  <  lO'*^  and  10  <  P  <  .999999,  procedures  were  given  for  computing  P  or  its  inverse  R 

(where  P  is  given  in  place  of  R)  to  6  decimal  -  digit  accuracy.  Here  procedures  are  given  for  computing 
P  or  R  to  6  significant  digits,  when  inherent  error  is  negligible,  for 

4-  4-  <  1/eps,  10  "  “  <  P  <  .9999999,  (8) 

where  eps  is  the  smallest  positive  machine  dependent  number  such  that  1  4-eps  >  1.  The  double  pre¬ 
cision  value  of  eps  for  the  IBM  PC  is  2  ~  ~  2.22  •  10  “  in  single  precision  eps=  2  ~  ~  1.19  10“^. 

Moreover,  for  the  smaller  ranges  used  in  [8],  P  or  R  can  now  be  found  to  approximately  8  significant 
digits.  All  computations  for  testing  were  carried  out  in  double  precision  PC  Fortran  on  a  Compaq 
Deskpro  386/20.  The  portable  double  precision  Fortran  function  which  yields  P,  given  R,  H,  K,  L,  u, 
v,  w,  is  called  ELLCOV;  the  portable  double  precision  subroutine  which  outputs  R,  given  P,  H,  K,  L  , 
u,  V,  w,  is  called  ELINV3.  It  is  anticipated  that  the  single  precision  version  of  these  subprograms  will 
be  included  in  the  Naval  Surf^lce  Warfare  Center  (NAVSWC)  Library  of  Mathematics  Subroutines  [13]. 
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n.  COMPUTATION  OF  P 

Initially  in  the  computation  of  P  by  ELLCOV  four  tests  are  used  to  eliminate  most  cases  where 
P  <  Z4  =  max(10  ~  100  epsm)  or  P  >  1  —  Ej,  where  epsm  is  the  smallest  machine  dependent  positive 

number  the  computer  can  use,  and  Ej  =max(10~**,  50  eps).  In  single  precision  for  the  IBM  PC, 
epsm  ~  1.17  •  10  ~  for  double  precision  epsm  ~  2.22  •  10  “ 


Test  #1 


P  set  to  0  if  R*  <  1.5^^  u  V  w  Z4. 


For  Test  #2,  starting  from  (1)  with 


D  =  M  =  max(u,  V,  w), 

V  =  E[x,/(>f2  u)]  E[y,/(>f2  v)]  E[z,/i^  w)]. 


we  have,  for  D  >  R, 


Hence, 


V  <  E[xi/(>[2  M)]  E[yi/(>r2  M)1  E[z,/(>f2  M)]  <  E[(D  -  R)/{>l2  M)]. 


P  ^  2x^uvW  M)]  lx r3  <  Z4. 


Test  #2  ; 


P  set  to  0  if  R^  E(7)  <  1.5>j^  u  vw  Z4 
7  =  (D-R)/(MM2). 


Generally,  when  none  of  the  above  tests  are  satisfied,  P  is  computed  by  the  numerical  integration  of 
(7).  However,  there  are  three  situations  where  P  can  be  evaluated  without  resorting  to  quadratures,  and 
a  fourth  where  P  is  given  by  P^..  We  shall  call  these  cases  A,  B,  C,  and  D. 
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CASE  A: 


For  small  R 


The  equivalent  of  (I)  is  to  take  the  normal  distribution  with  mean  at  (H,  K,  L)  and  to  place 
the  target  sphere  SP  at  the  origin.  Then  using  spherical  coordinates 

Xj  =  p  cos^  sinc^,  0<p<R,  0<^<  2r,  0  <  ^  <  t 

yj  =  p  sinO  sin^ 

Zj  =  p  cos^ 

with  the  volume  element  given  by  dV  =  p^  sin^  dp  d^  d^,  we  obtain 


P  = 


1 _ [ ^  f  f ""  E  cost?  sm_.^_-_H\  ^  /p  sing  sin<i  -  Kx  ^  /p  cc«^  -  L\ 

25r-^uvwjQjo  Jov  ^  J  \  ^{2  w  J 


Expanding  each  of  the  exponentials  about  p  =  0  and  carrying  out  the  integrations  gives,  after  some 
tedious  algebra, 


2-|-k2-(-l2) 


(9) 


where  H,  K,  L  and  Rj,  Rj'  defined  in  (5)  with 

T  =  2[Ri2(2H2-1)-1-R2^(2K2-1)-1-R32(2L2-1)] 

^r  =  2k{'^'~®[V(4H2-l)-hR2^4K2-l)-hR34(4L2-l)]}. 

Then  (9)  is  used  to  compute  P  when 

^{t2-(-8|  Ri‘‘(4H2-1)-(-R2'‘(4K2-1)-(-R3<(4L*-1)|  }<max  [5- 10’®,  eps]. 

Summarized  results  for  cases  B  and  C  follow.  The  derivations  of  the  final  expressions  either 
have  been  given  in  referenced  papers  or  they  will  be  given  in  the  Appendix  A  of  this  report. 


CASE  B: 


For  u  =  V  =  w  I 


P=iaerf(D,R)-;^R 


1  — e 


-4RDl 


4RD 


Je(D-R),  D  =  >|H*  +  k2-|-L2  ^0  (10) 

(11) 


P  =erfR-;^  R  E(R)  ,  D  =  0, 
where  E(z)  is  defined  in  (2),  R  and  L  are  defined  in  (5),  and 

fX 

aerf(  D,  R)  =  erf(D  R)  -erf(D  -  R),  erf  x  =  ^  E(t) 

An  algorithm  is  given  in  [9,  App.  A]  which  yields  aerf  to  13  significant  digits.  Equation(lO)  is  derived 
in  [8,  App.  A]. 


dt. 


(12) 
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Subtraction  of  leading  positive  terms  in  (10)  and  (11)  can  result  in  loss  of  accuracy  when 
computing  P.  Consider  (11)  for  small  R.  Note  that 

—  2  p/R\  V'  _ 2° _ R  2n  -I- 1 


erfR  =  AE(R)^-  l^R- 


[1,  p.297] 


so  that  the  leading  term  cancels  the  2nd  term  in  (11).  Therefore  when  R  <  .071,  (11)  is  replaced  by 


1 


3---(2n-(-l) 


R 


2n-t-l 


D  =  0. 


(13) 


In  the  case  of  (10)  there  are  two  situations  where  accuracy  may  be  lost.  First, 
consider  D  —  R  large  and  4 RD  >  —  log  eps  (For  the  IBM  PC,  log  eps  ~  —  36.044).  Then  with 

aerf(  D,  R)  =  erfc(D  —  R)  —  erfc(D  +  R),  erfc  x  =  1  —  erf  x,  (14) 

E(x) 


use 


erfc  X  ~ 


Vir  X 


1+ 

Y  (2x2)” 


(x-»oo),  [1,  p.298]  (15) 


in(14).  It  is  easy  to  see  that  the  term  erfc(D  -bR)  can  be  dropped  since  e  is  negligible  compared 

to  one.  Therefore  (10)  becomes 


-  2>fW 


or 


1  I  1  Vr  n"  •••  (2n-l) 

d-r'^d-rZ-^  ’  rom-Ri*!" 


p  ~ 


r  [2(d-r)^]' 

E(D-R)  l-3-(2n-l) 

D  4-^  ^  [2(D-R)2]" 


2>nr  (D-R) 


1 


(16) 


The  relation  (16)  for  P  is  used  when  D  —  R  >  4.25,  4RD  >  —  log  epw  ,  and  R  >  0.425  . 

Loss  of  digits  can  also  occur  when  R  is  small  and  D  >  0  in  much  the  same  way  as  was  seen  in 
obtaining  (13).  Here  we  use 

p_  4  p/p\y^  "  **2n-l(^)/'^  p2n  +  l 
P-;jfE(D)^  (2„  +  l)!  ^  ’ 


1 


(17) 


where  Hj^(x)  denotes  the  Hermite  polynomial  of  degree  k  [1,  p.771].  Equation  (17)  is  used  when 
R  <  .425.  It  is  derived  in  Appendix  A. 

Subprogram  ELLCOV  calls  subprogram  EQSIG  to  evaluate  P  from  the  above  relations  when  it 
recognizes  Case  B  from  its  input. 
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CASE  C: 

Let 

Z  =  ^|l-(w/u)2|  ,  S  =  L-RZ2,  F  =  L  +  RZ2. 

Using  the  fact,  which  is  easily  shown,  that 

Pp(f,0,0,u,u)  =  1 -E(rw/u),  r  =  -Jr^-(L-  z)^,  [5,6,7]  (18) 

and  substituting  this  result  into  (6)  gives  separate  relations  for  P  when  u  >  w  and  w  >  u: 

u_>^  (i  =  ^  —  1  ) 

P  =  iaerf(L,R)-^exp{-(f)2[R2-(L/Z)*]}aerf(|,RZ),  L  7^  0,  S<0  (19) 

P  =  iaerf(L,R)-^^^^[E(i|)erfc(|)-e-^«-I^E(i|)erfc(|)],  L  ^  0,  S>0  (20) 

P  =  erfR-RE(Rw/u)  L  =  0.  (21) 

w  >  u 

P  =  iaerf(L,R)-^^^^[daw(|)-e-4RLdaw(|)],  L^O  (22) 

P  =  erfR-:;|RE(R)^^^^|^,  L  =  0,  (23) 

where  ^ 

daw(x)  =  E(x)  J"^  e^  dt  (Dawson’s  integral)  .  [1,  p.298  ]  and  [8]  (24) 

Case  C  uses  (19)  — (23).  They  were  derived  in  [8,  pp.  A4  — A7].  In  a  number  of  situations 
however  these  equations  are  written  in  different  forms  here  to  reduce  the  computational  loss  in 
accuracy.  The  modified  forms  appear  below  in  (25)  — (31)  with  their  derivations  given  in  Appendix  A. 
In  each  of  the  modified  relations  given  below  the  conditions  under  which  they  are  used  are  spiecified 
first.  The  quantity  P  for  case  B  is  denoted  here  by  Pb-  and  Hn(x)  denotes  the  Hermite  polynomial  of 
degree  n.  These  results  are  used  in  the  subroutine  SEQHZ3  which  is  called  by  ELLCOV  when  the  latter 
has  recognized  case  C  from  the  input.  A  flowchart  for  SEQHZ3  is  given  at  the  end  of  this  section. 

LET:  {  4RL  <  10,  and  R  <  •^,  and  Rw/u  <  1. 

P=  2  R(R^/u)2E(L/>f2)^F2„^.,  (25) 

1 

^2n  +  1  —  2n  +  1  [^2"  - 1  ~  2  (Rw/u)^  Fj^  _  j]  ,  n  >  1 


Fi  =0,  Ho  =  1,  Gi  =  2  E(L/>[2)  . 
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LET: 


R  >  or  Rw/u  >  1,  and 
L  =  0  and  R  Z  <  '43. 

P  =  P  B  -  ^  E  ( R )  S  ( R  Z )  2" 


u  >  w 


w  >  u. 


(26) 

(27) 


LET: 


4RL  >10  or  R  >  ^  or  Rw/u  >  1,  and 
■  L^O  or  RZ>>l3,  and 

w/u>l/10  or  (w/u)^max(R, L)  >  1/2. 


Let 


erfcr(x)  =  1/^  —  X  erfc  X  ,  x>4 


efsz  =  F  erfcr(S/Z)  -S  E(2a|^)  erfcr(F/Z)  . 


Then 

E(R-L)rRz2 

2FS 

4  R^Z^  +  ( 1  +  E(  24^ ) 

—  efsz|, 

u  >  w,  1  >  5, 

(28) 

and 

^  =  ^B  -  ^4llL^^- -  R  [1  +  E(  24RL )]  -h  bdaw  1 }, 

w  >  u,  2 

(29) 

where 

l/(2y)  [1 +  Frac(y)/y^]  ,  y>5 

daw( y ) ~ 

[2] 

bdawl  =S/f2  Frac(F/Z)-F/S2  E(2>(RL)Frac(S/Z). 


The  computation  of  erfcr(x)  is  carried  out  by  the  Fortran  function  ERFCR  given  in  [13].  The 
evaluation  of  daw(x)  is  based  on  minimax  rational  approximations  [2].  For  t  >  5,  the  quantity  Frac(t) 
represents  a  rational  function  in  1/t^.  The  function  bdawl  is  computed  by  the  Fortran  subroutine 
BDAWl  given  in  Appendix  B. 


LET: 


Let 


4RL  >10  or  R  >  or  Rw/u  >  1,  and 
L  =  0  or  R  Z  >  •45,  and 
w/u  >  1/10  or  (w/u)*max(R,  L)  >  1/2. 


dxdaw(x)  =  daw(x)/x. 


Then 

P  =  1/2  aerf(L,R)-(2/Vlf)RE(L-R)dxdaw(RZ),  w>u,  L  =  0,  |  <  5, 

where  dxdaw(x)  is  computed  by  the  Fortran  function  DXDAW  which  is  given  in  Appendix  B. 


(30) 
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LET: 


4RL  >10  or  R  >  ^  or  Rw/u  >  1,  and 
<  L  ^  0  or  R  Z  >  and 

w/u  <  1/10  or  (w/u)^max(R,  L)  <  1/2  . 


With  i  =  -  1  , 

P  =  l|(l  -  8  )>ert(L,R) 

_ ^  - (w/„)“ (R^ - L  ^  Hn(-iL»/")  ^11'  +  ,4w,n  , n  j, 

The  recurrence  relations  used  to  generate  (31)  in  SEQHZ3  are  given  in  Appendix  A. 


(31) 


CASE  D: 


For  L  —  R  <  —  6  and  r  = 


r2-(L 


-A{2Wt) 


^~R 


P  =  Pc(R,H,  K,u,v). 


(32) 


This  case  is  recognized  when 

(  L  —  ■'(2  w  A  )^  <  R^min(eps/3,  10  ”  *^),  (33) 

where  A  replaces  the  lower  limit  of  integration  in  (7).  For  example,  let  R  =  10^°,  u  =  10'°,  v  =  2, 
w  =  1,  H  =  10'°,  K  =  1,  L  =  2.  In  this  case  L  —  R  =  (2  —  10 *°)/>j2  can  safely  be  replaced  by 
A  =  -  7.0  since  r  in  (7)  is  an  increasing  function  of  t.  Then 


(L-V2wA)VR^=;1210"^°<  (2.22/3)  •  10  “  and  erf(7.0)  ~  1  -  4  •  10  ' 


Hence  P  =  P^.(  R,  H,  K,  u,  v )  =  .47724  986805. 


In  general,  excluding  cases  A,B,C,D,  the  probability  P  is  computed  from  (7)  by  numerical 
integration.  A  24th  order  Gaussian  quadrature  rule  [I,  p.916]  is  used  for  this  purpose.  The  primary 
objectives  are:  (a)  to  determine  the  order  of  integration,  i.e.,  which  integration  should  be  carried  out 
last  as  the  t-integration  in  (7),  (b)  to  obtain  effective  limits  of  integration. 


From.  Test#3,  T  <  9.60.  Indeed,  if  L  —  R  >  9.6,  then 


(.  oo  _  oo 

P<:^l  E(z)  Pc(r,H,K,u,v)  dz  <:i  E(z)  dz  =  ierfc(9.6)  ~  2.8  •  10 
^  9.6  ^  9.6 

In  addition,  the  lower  limit  in  (6)  or  (7)  can  be  restricted  to  values  greater  than  —  7.0.  In  fact,  by 
Test#4  and  (10),  when  R>D  ,  then  P  >  P*  ~  aerf(D,  R)/2  where  the  second  term  in  (10)  can  be 
neglected  for  P  —  D  >  7'j2  M.  Therefore  since  P  is  no  less  than  P*, 

P  >  aerf(D,R)/2  >  erf(7.0)  ~  1  -  4  ■  10  “ 
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Consequently,  if  we  identify  an  initial  effective  integration  range  in  (7)  to  be  from  A  to  B,  then 

A  =  max(L  -  R,  -  7.0),  B  =  min(L,  9.6).  (34) 

It  is  understood  that  an  A  and  B  are  determined  for  the  three  different  integration  orders,  i.e.,  with  the 
original  x,y,  and  the  z-integrations  reordered  so  that  each  is  carried  out  last. 

Let  the  integrand  of  (7)  be  denoted  by  G(t).  Then  ELLCOV  calls  TQUAl  to  determine 
whether  A  can  be  raised  and/or  B  lowered  from  the  initial  values  given  in  (34).  At  each  Gaussian 
abscissa  tj  on  [A,B]  the  integrand  of  (7),  G(tj),  is  evaluated  at  increasing  i  starting  at  i  =  1,  where 
A  =  t(j  <  tj  <  .  .  .  <  t23  <  t24  <  tjs  =  B.  Set  the  quantity  T8  =  max  (100  epsm,  10  ~  ^^),  with  epsm  as 
defined  on  page  3.  Then  for  the  smallest  t^  for  which  G(tj)  >  T8,  A  is  replaced  by  t,  _j.  Similarly, 
starting  at  j  =  24  and  with  decreasing  j,  the  largest  t^-  for  which  G(tj)  >  T8  is  found.  B  is  then  replaced 
t>y  +  function  G(t)  is  computed  by  the  Fortran  function  FN2  given  in  Appendix  B. 

At  this  point  three  sets  of  integration  limits  A  and  B  have  been  determined,  one  for  each 
different  l2ist  integration.  Indicate  these  limits  by  AX,  BX;  AY,  BY;  AZ,  BZ.  The  final  order  of  inte¬ 
gration  is  now  selected  as  the  one  for  which  the  integration  interval  (B  —  A)  is  the  largest.  An  exception 
to  this  choice  is  made  if  one  lower  limit  of  integration  is  <  —  2  and  the  other  two  lower  limits  are  not 
negative,  for  example,  AX  <  -  2  and  AY  and  AZ  are  not  negative.  In  this  particular  case  the  x 
integration  is  chosen  as  the  last  integration.  If  the  above  exception  does  not  hold  and  a  tie  occurs,  say 
between  the  x  and  y  integrations,  (BX  —  AX  =  BY  -  AY)  the  x-integration  is  carried  out  last  if 
AX  >  AY  and  AX  <  —  2,  otherwise  the  y-integration  is  carried  out  last.  Choosing  the  order  of 
integration  in  the  ways  described  above  is  based  on  numerical  studies  and  the  heuristic  argument  that 
the  more  spike-like  the  integrand  the  more  difficult  it  is  to  obtain  an  accurate  numerical  integration. 

With  the  appropriate  interchanges  having  been  made,  the  integration  to  be  carried  out  last  has 
now  been  set  and  it  is  indicated  by  the  z  or  t-integration  as  shown  in  (6)  or  (7),  .  Two  further  attempts 
are  subsequently  made  to  improve  the  values  of  A  and  B. 

It  may  occur  that  for  some  t  in  [A,B],  ~  1  —  c,  t  =  2.3  •  10  “  in  which  case,  since  r  in  (7) 

is  an  increasing  function  of  t, 

P  ~  I  [  E(t)  -1-  E(2  L  - 1)]  P J  r,  H,  K,  u,  v )  d  t  -»-  i aerf(L,r), 

A 

where  Y  is  defined  below.  Such  a  situation  is  recognized  by  employing  previous  results  fron.  [3,  p.l5] 
and  [4,7].  We  have  P^  >  1  -  2.3  •  10  “  for  t  >  t,  if 

r  =  ,  R^  -  (L  -  V2  w  t)  ^  >  ■>  H  ^ -f- K  ^ -f  7m  =  G,  m  =  max(u,v). 
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or 

|L->[2wt|  <a1r2-  =r, 

or 

t  =  L-r<i<L  +  r,  r  =  r/(>|2w). 

At  this  point  a  final  effort  is  made  to  improve  A  and  B.  ELLCOV  calls  subroutine  SQUAD 
which  first  evaluates  the  integrand  in  (7),  G(t),  at  each  Gaussian  abscissa  on  [A,  B]  starting  from  the 
endpwints  and  moving  symmetrically  toward  the  center  of  (A,  B].  An  estimate  is  now  available  for  P. 
With  this  estimate,  the  same  procedure  used  in  TQUAl  is  carried  out  by  SQUAD  to  jwssibly  further 
improve  A  and/or  B.  The  difference  here  is  that  the  better  estimate  for  P  replaces  T8  used  in  TQUAl. 
If  however  both  A  and  B  are  unchanged  by  SQUAD  then  the  value  for  P  obtained  from  SQUAD  gives 
the  final  result  for  P. 

If  A  and/or  B  are  improved  by  SQUAD,  then  ELLCOV  calls  subroutine  RQUAD  to  obtain  P 
by  a  final  24th  order  Gaussian  numerical  integration,  based  on  the  latest  values  of  A  and  B  found  from 
SQUAD.  In  SQUAD  or  RQUAD,  we  have 


where 


p  =  ^^Ey(*)[G(t7)+G(tt)], 

i  =  12 

G(t)  =  [E(t)  E(2L  - 1)]  Pc[T(t),H,K,u,v] 


r(t)  =  ,  R^  -  (L-'^wt)^ 


‘7 


B-l- A 

2 


-x(i) 


B-A 

2 


B-l- A 


2 


+  x(i) 


B-A 

2 


y(i)  =  24th  order  Gaussian  weights  on  [  —  1, 1] 

x(i)  =  24th  order  Gaussian  abscissae  on  [-  1, 1].  [1,  p.  916] 

The  x(i)  and  y(i)  are  stored  in  ELLCOV  (see  listing  in  Appendix  B). 

A  parameter  J  is  set  in  ELLCOV  to  either  0  or  1  or  2.  It  is  set  to  0  if,  referring  to  P^.  above, 
H  =  K  =  0;  it  is  set  to  I  if  u  =  v.  In  either  of  these  two  cases  the  subroutine  CIRCV  is  used  to 
evaluate  P^.,  otherwise  the  subroutine  T^'KILL  is  used  for  this  purpose  with  J  set  to  2.  It  is  advantageous 
to  use  CIRCV  if  possible  since  it  evaluates  P^  more  accurately  than  PKILL  and  is  about  ten  times 
faster.  However  the  order  of  the  integration  is  determined,  as  described  above,  by  TQUAl  and  no  effort 
is  made  to  modify  that  choice. 
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The  quantity  f(t,)  can  be  simplified  when  A  =  L  —  R  and  =  L.  In  this  case,  it  is  easy  to  show 

that 

f(t-)  =  RU(i)  =  R^l-[l-|-x(i)]V4  ,  i  =  12,ll . 2,1 

r(t+)  =  R  U(12  -(-  i)  =  R^l-[l-x(i)]V4  ,  i  =  12, 11,. .  .,2, 1. 

The  array  U(j)  has  been  precomputed  and  is  available,  thus  eliminating  the  computation  of  24  square 
roots  and  a  possible  loss  of  accuracy.  When  this  feature  comes  into  play,  the  parameter  IL  in  SQUAD 

is  set  to  1.  The  array  is  stored  in  FN2  where  G(t)  is  computed,  (See  the  listing  in  Appendix  B). 


Although  we  claim  6  significant-digit  accuracy  for  P,  it  may  not  be  possible  to  obtain  this  ac¬ 
curacy  for  extreme  values  of  the  input  where  inherent  error  plays  a  dominant  role.  For  example,  using 


ELLCOV  on  an  IBM  PC  with  double  precision  Fortran  (  ~  16  significant  digits),  let 
H  =  10l^  K=0,  L  =  0,  u  =  5,  v  =  w=  1, 

then 


R 

P 

1.0000  00000  0026-10*^ 

.99999  99004 

1.0000  00000  0025  10*^ 

.99999  97136 

1.0000  00000  0000-10^^ 

.50000  00000 

9.9999  99999  9999  •  lO'^ 

.49214  92357 

9.9999  99999  9600  •  lO'* 

.62292  51932-10-^® 

9.9999  99999  9599  •  lo'^ 

.52932  23990  10"^^, 

The  underlined  digits  are  in  doubt. 

A  few  more  numerical  results  are  included  where  inherent  error  is  not  a  factor.  If 
H  =  1.0,  K  =  2.0,  L  =  3.0,  u  =  2,  v  =  4,  w  =  1, 

then 

R  P 


1.0000  00000  0000- 10“^ 

.28764  95875- 10"® 

5.0000  00000  0000  - 10 

.43140  77476  - 10”^ 

1.0000  00000  0000 

.53583  15121-10“^ 

3.0000  00000  0000 

.79694  49445  - 10  ~  ’ 

4.0000  00000  0000 

.25115  76137 

5.0000  00000  0000 

.47295  03432 

7.0000  00000  0000 

.78979  21674 

1.2000  00000  0000  10^ 

.98945  55844 

1.8000  00000  0000  10’ 

.99994  81467  . 
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We  also  take  the  opportunity  here  to  correct  3  typographical  errors  in  Table  1  of  [8,  p.l2]. 


R=  1 

H  =  1/2 

K  =  L  =  0 

c 

II 

< 

II 

II 

P  =  .  17955  97978 

R  =  1 

H  =  2 

K  =  L  =  0 

u  =  V  =  w  =  ^2/3 

P-. 33473  93607-  10-' 

R=  1 

H  =  0 

II 

c-l 

II 

o 

u  =  2v  =  2w  =  2/^j3 

P  =  .37539  30077 

Extensive  testing  was  carried  out  to  establish  the  accuracy  claimed  above  for  the  computation 
of  P  by  ELLCOV.  This  testing,  using  a  Compaq  Deskpro  386/20  PC,  compared  ELLCOV  double  pre¬ 
cision  results,  with  the  results  from  an  adaptive  quadrature  subroutine,  DQAGS,  which  is  contained  in 
NSWCLIB  [13]. 
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FLOWCHART  FOR  SEQHZ3 

CASE C:  H  =  K  =  0.  u  =  v w.  Ws  (w/u)^ .  Y  s  1 1  -  W |  ,  S  h  L -RY 
/  START  y 


4RL>  10  or  r2>2  or  R^  >  1? 


|yes 

pH 

InO 

V  CO 

I  >  11/ 

4  “ 

r— 

L)t0  or 

W<.01  and 

V 

R2y  >  3? 

max  (L,  R)  W<l/2? 

15 


15 


■  u^  w? 

'YES 

L  = 

0? 

1 

KS 

SSO?  - 

^YES 


•  Pn  =  P  for  u  =  V  =  w. 

D 

•  E  =  max  (eps  /  2, 10  )  See  page  2  for  eps. 

•  KKK:  integer  used  to  identify  a  particular  branch  of  SEQHZ3. 

•  Numbers  in  rectangles  refer  to  equation  numbers  of  the  text 

•  Numbers  outside  boxes  and  circles  refer  to  Fortran  labels  of  the  source  code. 
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ra.  COMPUTATION  OF  R  (The  inverse  problem) 

In  this  section  a  procedure  is  described,  given  P,  H,  K,  L,  u,  v,  w,  to  estimate  R.  The  method 
for  estimating  R  within  a  given  accuracy  generally  requires  a  sequence  of  values  of  P  from  ELLCOV. 
Since  these  computations  usually  involve  time-consuming  numerical  triple  integrations  an  effort  is 
made  to  keep  them  to  a  minimum  by  obtaining  a  good  early  estimate  of  R  .  Once  such  an  estimate  is 
found,  halving,  regula-falsi,  and  King’s  root  finding  procedure  [12]  are  used  to  obtain  a  final  estimate, 
R  ,  for  R  .  The  objective  is  to  obtain  R  to  at  least  six  significant  digits  (when  inherent  error  does  not 
contaminate  the  computations)  for  the  range  specified  earlier  in  (8),  namely 

H^-t-K^-l-L^  <  1/eps,  10"“  <  P  <  .9999999.  (35) 

The  Fortran  subroutine  ELINV3  yields  R  .  A  listing  of  the  code  is  given  in  Appendix  C  with 
listings  for  three  additional  portable  supporting  subprograms  ELLRC,  FCNl,  SUB3. 


The  unknown  value  of  R  will  always  be  contained  in  a  known  interval  [R^^,  Rf/]-  Initial  values 
of  the  lower  and  upper  bounds,  R^^  and  R;^,  are  found  for  R  . 

Let  I  =  max  [3.Jx/2  Puvw,  (,Jir/2  P  M)^]  ,  M  =  max(u,v,w).  Then 


D,  ifP>l/2andI<D^  D  = 
otherwise. 


(36) 


The  derivation  of  (36)  is  discussed  in  [8].  An  initial  value  of  R^  is  obtained  from  the  relation 
R<  R// =  D-(-B(j)M  with  A(j)  =  P(B(j)M,0,0,0,M,M,M), 
where  B(j)M  gives  the  radius  of  a  sphere  centered  at  the  origin  such  that  A(j)  is  the  smallest  value  for 
which  A(j)  >  P,  This  approximation  is  discussed  more  fully  in  [8;  pp.  14-15].  The  arrays  A(k)  and  B(k) 


are  given  by 


A(l)  = 

10 

-30 

B(l) 

=  1.56- 

10- 

-  10 

A(2)  = 

10 

-  25 

B(2) 

=  7.23  • 

10- 

-  9 

A(3)  = 

10 

-  20 

B(3) 

=  3.36  • 

10' 

-  7 

A(4)  = 

10 

—  15 

B(4) 

=  1.56- 

10' 

-  5 

A(5)  = 

10 

-  10 

B(5) 

=  7.22  • 

10' 

-4 

A(6)  = 

10 

-  8 

B(6) 

=  3.36  • 

10- 

-3 

A(7)  = 

5- 

10-® 

B(7) 

=  2.66  • 

10- 

-  2 

A(8)  = 

10 

-  4 

B(8) 

=  7.23  • 

10- 

-  2 

A(9)  = 

10 

-  2 

B(9) 

=  0.339 

A(I0)  =  0.10  B(10)  =  0.765 

A(ll)  =  0.30  B(ll)=  1.1933 

A(12)  =  0.60  B(12)  =  1.7170 

A(13)  =  0.90  B(13)  =  2.5005 

A(14)  =  0.999  B(14)  =  4.0335 

A(  15)  =  0.999999  B(15)  =  5.5380 

A(  16)  =  0.99999999  B(16)  =  6.3500 

A(17)  =  1.00  B(17)  =  7.7000 


14 


NAVSWC  TR  91-487 


An  improvement  for  R  is  generally  obtained  by  using  an  estimate  R^  given  by  F.E.  Grubbs 
[11].  His  estimate  depends  on  the  percentage  point  of  the  Chi-squared  distribution  which  is  determined 
by  using  the  subroutine  GAMINV  contained  in  the  NAVSWC  math  library  [13].  The  quantity  actually 
given  by  Grubbs  is  R(j^,  namely 

Rg'  =  X(4/V5,  P)  -  4/V5  ]  W2  -h  Z},  (37) 

where 

=  u^  -1-  v^  -I-  w^ 

Z  =  l-^-DV<r^ 

V4  =  2{i4[H-4H2]-h3il[n.4K2]+  h1[i+4l2]| 

^  (T  u  (T  ^ 

T5  =  8  {  4  [1  +6H2]  +  6K']  +  ^  [1  +  6L2]  } 

^  (T  (T  (T  ' 

W2  =  T5/(2V4),  Vs^TsVv,^ 
and  X(A,P)  satisfies  the  incomplete  gamma  function  relation 

a  =  4/V5.  [l,p.255] 

The  quantity  R^j^  from  (37)  occasionally  may  give  a  poor  estimate  for  R*  or  it  may  even  give  a 
negative  value  and  consequently  be  of  no  use  in  such  cases. 

We  attempt  to  improve  our  initial  estimate  for  R  by  finding  a  constant  R(;7,  where 

r  L  -f  R^  ^  K  -i"  R^  p  H  -(-  R^ 

P=  __  __  __  &’(xi,yi,Zi,u,v,w)  dxj  dyj  dzj  .  (38) 

J  L  -  R^  K  -  Rc  ■*  H  -  Rc 

The  quantity  S'  is  defined  in  (2).  The  right  hand  side  of  (38)  gives  the  cumulative  normal  probability 
P  over  the  cube,  with  center  at  (H,  K,  L)  and  equal  sides  parallel  to  the  Xj,  yj,  and  Zj  axes,  which 
contains  the  sphere  with  the  same  center,  and  of  radius  R^^.  The  subroutine  ELLRC,  using  Newton- 
Raphson,  gives  R^.  Hence  R^  <  R  <  ^3  R^.  The  maximum  number  of  Newton-Raphson  iterations  al¬ 
lowed  by  ELLRC  is  40.  The  Fortran  listing  for  ELLRC  is  given  in  Appendix  C. 

With  the  latest  values  for  R^  and  R^  and  the  corresponding  P  values  P^^  and  P//,  a  halving 
procedure  is  carried  out  until  Pjl>10“^  P^  and  (R;^  —  R^)  <  R/10.  When  both  of  the  above 
inequalities  are  satisfied  ELINV3  proceeds  to  obtain  a  final  estimate  for  R  by  employing  King’s  method 
[8,  12].  The  method  was  described  with  a  flow  chart  in  [8].  It  may  be  looked  upon  as  a  modified  regula- 
falsi  procedure.  An  estimate  R  for  R  is  considered  satisfactory  if  |P(R  )  -  Pj  <  E8,  where  E8  depends 
on  P  and  is  given  by 
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max(10  eps,  10  *)  P  0  <  P  <  .999 

max(10  eps,  10”^)  .999  <  P  <  .999999 

E8  =  <  (See  page  2  for  eps). 

max(10  eps,  10  " .999999  <  P  <  .99999999 

max(10  eps,  10  ~  .99999999  <  P  <  1.0. 

If  the  above  inequality  |P(R  )  —  P|  <  E8  is  satisfied  an  indicator  IND  is  set  to  0  in  ELINV3.  Other  exits 
from  ELINV3  are  also  identified  by  IND,  namely, 

IND  =  -  1  P  <  max(10  10®  epsm)  (See  page  3  for  epsm). 

IND=-(-l  P>  l-max(10"*^  lOeps) 

IND  =  -1-2  N>30 

IND  =  -1-3  R;^  -  R^^  <  max(10  eps,  10“^^)R 

IND  =  -f  4  ELLCOV  not  able  to  yield  P(R  )  with  the  accuracy  required. 

If  IND  =  |1|,  P  has  been  specified  outside  allowable  limits.  The  output  R  is  set  to  —  10^®.  Let  N  denote 
the  number  of  iterations.  If  IND  =  2,  then  the  maximum  allowable  number  of  iterations  (30)  for 
finding  R  was  reached.  N  is  increased  by  one  for  each  call  to  ELLCOV  for  which  its  output  is  nonzero. 
Extensive  checking  has  found  N  <  27.  If  IND  =  3,  then  more  digits  are  required  in  the  R  estimates 
than  are  available,  .i.e,  the  word  length  for  the  particular  machine  in  use  is  too  short.  If  IND  =  4,  then 
the  subprogram  that  yields  P(  R  ),  ELLCOV,  cannot  produce  the  accuracy  required.  This  difficulty  can 
only  appear  if  P  is  very  near  one  and  it  is  due  to  either  the  limitations  of  the  numerical  integration  or 
to  inherent  error. 

Sonie  numerical  results  from  ELINV3  are  given  below  in  Table  1  (4  parts)  which  reproduces 
Table  2  of  [8]  with  improved  estimates  for  R .  The  notation  for  Table  1  follows.  Let  P  denote  the  input 
values  of  P  and  let  R  denote  the  final  estimate  for  R  using  ELINV3.  RERR  represents  the  relative  error 
given  by  1 1'(  R  )  —  P  [/  P  .  IND  and  N  have  been  defined  above. 
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TABLE  1  (PART  1) 


u  =  2.0 

V  =  4.0 

w  = 

6.0 

H  =  5.00 

K  =  10.0 

L  = 

20.0 

R 

RERR 

IND 

N 

P 

0.3183274539-10' 

-(-0.260-10"* 

0 

8 

0.500000  - 10  ■  ® 

0.1028906156 -10* 

-  0.392  - 10  ~ 

0 

8 

0.500000-10"* 

0.1660907596  10^ 

-0.641-10"'* 

0 

6 

0.100000 

0.2339931591  •  10* 

-(-0.163-10"® 

0 

5 

0.500000 

0.3048621773-10* 

-(-0.708-10"* 

0 

7 

0.900000 

0.4077123752-10* 

-0.231-10"'® 

0 

7 

0.999000 

0.4844201564-10* 

-(-0.956-10"'® 

0 

6 

0.999995 

TABLE  1  (PART  2) 


u  =  5.0-10"® 

V  =  5.0-10" 

2 

w  = 

1.0 

H  =  0.20 

K  =  1.00 

L  = 

5.00 

R 

RERR 

IND 

N 

P 

0.1171110995-10' 

-(-0.653 -10"® 

0 

9 

0.590000-10"® 

0.2629994892  -  lO' 

-0.127-10"'* 

0 

8 

0.500000-10"* 

0.3855992354  -  lO' 

-0.143-10"" 

0 

6 

0.100000 

0.5103195068-10' 

-(-0.287-10"'* 

0 

7 

0.500000 

0.6364036803  -  lO' 

+  0.572-10"* 

0 

7 

0.900000 

0.8154468511-10' 

-0.259-10"'® 

0 

7 

0.999000 

0.9472449895  -  lO' 

+  0.919-10"® 

0 

7 

0.999995 
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TABLE  1  (PART  3) 


u  =  2.0 

v  =  5.0 

W  = 

1.0 

If  =  0.00 

K  =  5.00 

L  = 

10.0 

R 

RERR 

IND 

N 

P 

0.6184610273  10^ 

-0.780-10"^ 

0 

10 

0.500000  - 10  "  * 

0.8185500434-10^ 

1 

p 

o 

o 

o 

o 

0 

13 

0.500000-10"^ 

0.9733490650  •  10^ 

-0.119-10-'^ 

0 

6 

0.100000 

0.1172659410-10^ 

+  0.211-10"'° 

0 

6 

0.500000 

0.1545659369- 10^ 

+  0.910-10"* 

0 

6 

0.900000 

0.2295500292  - 10^ 

+  0.285-10"'^ 

0 

8 

0.999000 

0.2902552953  - 10^ 

+  0.241-10"° 

0 

6 

0.999995 

TABLE  1  (PART  4) 


u  =  2.0-10'* 

V  =  2.0  - 10° 

w  = 

1.0 

H  =  0.0 

K  =  2.0  - 10® 

L  = 

0.0 

R 

RERR 

IND 

N 

P 

0.2137608898-10® 

-0.379-10"® 

0 

6 

0.500000-10"® 

0.1253323942-10* 

-0.374-10"* 

0 

1 

0.500000 -10  "2 

0.2513226958-10° 

-0.580-10"° 

0 

1 

0.100000 

0.1348979507-10*° 

+  0.942-10"*° 

0 

1 

0.500000 

0.3289707270  - 10*° 

+  0.697-10"** 

0 

1 

0.900000 

0.6581053496-10*° 

+  0.541-10"*° 

0 

1 

0.999000 

0.9129575506-10*° 

+  0.222-10"*® 

0 

1 

0.999995 
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IV.  TEST  RESULTS  FOR  ELLCOV  AND  ELINV3 

Extensive  testing  was  carried  out  to  establish  the  accuracy  claimed  for  ELLCOV  and  ELINV3 
keeping  in  mind  that  inherent  error  plays  a  major  role  when  D  is  large  (see  page  11).  The  tests 
compared  ELLCOV  results  with  the  corresponding  outputs  from  a  double  precision  adaptive  quad¬ 
rature  subroutine,  DQAGS,  which  is  contained  in  NSWCLIB  [13].  For  ELINV3,  values  of  P  ranging 
from  10“^®  to  .9999999  were  given  with  extensive  ranges  of  the  variables  H,  K,  L,  u,  v  (w  =  1).  An 
estimate  R  for  the  true  value  R  was  found  using  ELIN V3.  ELLCOV  was  then  used  with  R  to  compute 
P  for  comparison  with  the  initial  given  value  of  P. 

Testing  using  a  Compaq  Deskpro  386/20  computer  with  a  CYRIX  numeric  coprocessor  showed 
that  the  average  time  to  compute  a  value  of  P  using  ELLCOV  was  ~  0.4  seconds.  The  average  time 
using  ELINV3  to  find  a  value  R  was  ~  1.1  seconds  requiring  on  the  average  approximately  6  iterations. 
The  maximum  observed  number  of  iterations,  27,  occurred  for  isolated  cases  in  the  range 
10  -  20  ^  p  £  10  ~  3  maximum  observed  time  for  EL1NV3  was  25  seconds,  also  occurring  for  small 
values  of  P.  All  testing  was  done  using  double  precision  Fortran  for  the  IBM  PC. 
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DERIVATION  OF  (17)  AND  (25) -(31) 


DERIVATION  OF  (17): 


Equation  (17)  holds  when  u  =  v  =  w  and  D  ^  0  and  is  given  by 
P  =  ^E(D)^  D  = 

where  Hjj(x)  denotes  the  classical  Hermite  polynomial  of  degree  k.  From  (10) 

P  =  ^erf(D  R)  -  erf(D  -  R)  -  ^  [E(D  -  R)  -  E(D  R)] }. 

From  [9,  p.  A-3] 

OO  11  /T)\  p2n  -f  1 

aerf(D,  R)  =  erf(D  +  R)  -  erf(D  -  R)  =  ^  E(D)  ^  ■ 

Using  the  property  of  the  Hermite  polynomials 

E(x)  =  (-!)”  E(x)  H„(x),  [1,  22.11.7] 

and  expanding  E(D  -t-  R)  about  R  =  0  gives 

OO  u  pn 

E(D-|-R)  =  E(D)53(-l)”-^i5-j^/ - . 

T\  ~  0 

With  these  results  A(2)  becomes 

„_lf  4  2  E(D)^H2„+,(D)R="  +  '1 

(2„  +  l)!  Mr  D  (2n  +  l)!  /’ 


P  =  ^  T  £  (1?^  |2  D  H,„(D)  -  H,„  ^  ,(D)] . 

n  =  o'*  ^ 

Making  use  of  the  recurrence  relation  for  Hermite  polynomials,  namely, 

H„  +  i(x)  =  2xH„(x)-2nH„_j(x),  H^  =  1,  Hj  =  2x,  [1,22.7.13] 

in  A(7)  yields  A(l). 


A(l) 

A(2) 

A(3) 

A(4) 

A(5) 

A(6) 

A(7) 

A(8) 


DERIVATION  OF  (25): 


Equation  (25)  is  given  by 

P=  2  f^(f^^/y)2E(L/>l2)^F2„  +  i  A(9) 

1 

-  1  =  (2n-m2n-2)!  ^ 

F,  =0,  Ho=l,  Gj  =  2  E(L/^l2)  .  A(9.3) 
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Let 


/■  L  -f  R  - 

f(R)=  E(z){l-< 

J  L-R  ^ 


-(w/u)"[R^-(L-z)''] 


dz. 


A(10) 


Then,  for  H  =  K  =  0  and  u  =  v,  we  have  from  (18)  and  (6) 

P  =  if(R).  A(ll) 

Our  aim  is  to  obtain  the  Taylor  series  for  f(R)  about  R  =  0  and  also  the  recurrence  relation  for 
computing  the  successive  terms  of  the  series.  From  A(10)  and  (12) 


f(R)=^aerf(L,R)-  [^^^E(z)e  dz. 

^  J  L-R 

Note  that  f(R)  is  an  odd  function  of  R  so  that  f^^"^(0)  =  0,  n  =  0,1,2,...  ,  where 


A(12) 


Thus 


Differentiating  f(R)  in  A(12)  gives 


f(")(R)=^f(R). 

^  f(2n  +  l)(o) 


A(13) 


i(^)2f(*)(R)  =  R[^^*^E(z)e  (L  z)^] 

^  L-R 

=  R[Q(R)-f(R)],  f(‘)(0)  =  0,  A(14) 

where,  using  A(3)  with  D  =  L,  we  have  introduced  the  notation 

OO  II  IT  \  pZn  +  1 

Q(R)  =  ^  aerf(L,  R)  =  2  E(L)  ^  [SeeA(3)].  A(15) 

Successive  differentiations  beginning  with  A(14)  give 

i  (^)2  f(2)(R)  =  Q(R)  _  f(R)  +  R  [q(‘)(R)  -  f(»(R)] 

i  (^)2  f(3)(R)  =  2  [Q<')(R)  -  f<'^(R)]  -(-  R  [Q^^\R)  -  f<^\R)] 
i  (^)2  f(^\R)  =  3  [q'^VR)  - f^"^(R)]  +  R  [Q‘^^(R)  -  f'^VR)] 


i  f<")(R)  =  (n-  1)  [Q<"-^)(R)-f<""*>(R)]-l-R  [Q<"~*\R)-f<""‘’(R)]. 

Hence,  it  follows  that 

1  (U  )2  f (2n  +  ^  2n  [Q<2"  -  *\o)  -  f "  *\o)] ,  A(  16) 
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where,  from  A(15), 

q(2"  -  i)(0)  =  2  E(L)  _2(L),  q(')(0)  -  2  E(L). 

Now  define  F2„  +  i  by 

2  R  (Rf )  2  E(L>f2)  F2„  + 1  =  R "  =  1^2,.  -  •  • 

Then,  from  A(ll),  A(13),  and  A(16) 

oo 

P  =  ^  R  (R "  E(L/>|2)  E  R2n  + 1  '  Fi  =  O' 

n  =  1 

where,  using  A(15)  and  A(16), 


P2n  + 1  -  2n 


and 


(2n-  1)!  (2n-l)! 

=  2TtT[^2n-i-2(Rf)'F2„_,],  n>l, 
2E(L/aI2)  r2"-2 


R 


2n-2 


2n  (2n  + 1) 


’2n  -  1  - 


H2„-2(L). 


Also 


(2n-l)! 

F,  =0,  Ho(L)  =  l,  Gi  =  2E(L/>f2), 
where  the  first  relation  follows  from  setting  n  =  0  in  A(16)  and  using  A(14). 


DERIVATION  OF  (26)  -  (27): 

When  H  =  K  =  L  =  0  and  u  =  v,  one  obtains  from  (21)  for  u  >  w 

P  =  erf  R  -  RE(R  w/u)  Z  =  ^|l-(w/u)2|  . 


I  hen  using 


in  A(18)  yields  (26), 


erf(RZ)=^E(RZ)E  y: 


RZ 

2  ** _  CO  2n  1 


3-"(2n-l- 1) 


(RZ)- 


P  = 


where 


Pg  =  erfR-^  R  E(R)  [See  (11)]. 
For  w  >  u,  one  obtains  from  (23) 


P  =erfR-=^  E(  R) 


2  R  wDx  daw(RZ) 


,  L  =  0  , 


where 


daw(x)  =  E(x)  e^  dt 


J  0 


[1,  p.297] 

u  >  w. 


MW  RZ 

(Dawson’s  integral).  [1,  p.298  ]  and  [8] 


A(16) 

A(17) 

A(17.1) 

A(17.2) 

A(17.3) 

A(18) 

A(19) 
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Also 


daw(x)  ^^dt=x[  E(xco80)  cos^  AO,  t  =  x  si 

J  0  •'o 

1  + 1  r  ’’’/^ 


sin^ 


=  V(_l)nxi2^f  eos2'*  +  led0 

0  J  0 


_  i\n  x^"  +  ^  2“  n! 

-  id  1.3.5---2n-H  • 

Using  this  result,  with  x  =  RZ,  in  A(19)  above  gives  (27), 


P  =  PB-^E(R)^,,3<-|>;^l3(Rzd" 


w  >  u. 


DERIVATION  OF  (28)  -  (29): 


For  H  =  K  =  0  and  u  =  v,  with  S  >  5Z,  we  derive  the  following  relations: 


P  =  Pn- 


E(R-L) 


2 


erfcrx  =  I/Vjt  — X  e^  erfc  x  ,  x>4 

efsz  =  F  erfcr(S/Z)  —  S  e  “  erfcr(F/Z), 


and 


4RL 


72  1 

+  ^  bdawl  w  >  u 


daw(y)  ~  l/(2y)  [1 +  Frac(y)/y^]  ,  y>5  [2] 

bdawl  =  S/f2  Frac(F/Z)-F/S^  E(2>fRL)Frac(S/Z), 


with 


Z  =  ^|l-(w/u)2|  ,  S  =  L-RZ^  F  =  L-t-RZ2. 

For  the  derivation  of  (28),  which  is  A(20),  first  consider  (20),  with  u  >  w, 

P=iaerf(L,R)---^^[E(i|)erfc(|)-e-'‘*^*^E(i|)erfc(|)],  L  ^  0,  S>C 

The  quantity  in  square  brackets  in  A(23)  can  now  be  written,  using  A(20.1)  and  A(20.2),  as 

Z[F-Se-^RLj  ^ 

>fWFS  Fscisz. 

Substituting  this  result  into  A(23)  and  also  adding  and  subtracting 


A-4 


A(20) 

A(20.1) 

A(20.2) 

A(21) 

A(21.1) 

A(21.2) 

A(22) 

A(23) 
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2  2(L2-R2z^) 

4RL  2FS 


A(24) 


gives,  after  using  A(22), 


p-P  +RSri+e 

2FS  Vjr4RL  ^ 


^  -  efsz|  + 


_  2(L^-R^Z^) 

W4RL  2FS 


where  Pg  is  given  by  (10).  Carrying  out  the  obvious  cancellation  in  A(25),  followed  by  some  minor 
rearrangements  gives  A(20). 


In  order  to  derive  (29)  above  begin  with  (22),  where  w  >  u,  .i.e., 
D_l„„r/rox  E(L-R)[■^_,F^  ^-4RLj^^,S 


P  =  i  aerf (L,  R)  -  [daw (|)  - 


(§)],  L^i 


The  quantity  in  square  brackets  in  A(26)  can  now  be  written,  using  A(21.1)  and  A(21.2),  as 
^|s-Fe-4RL^Z2|’^  Frac(F/Z)  -  e  “  ^  ^  ^  ^  Frac(S/Z)j  |. 
Substituting  this  result  into  A(26)  and  adding  and  subtracting  A(24)  gives,  using  A(22), 


E(L-R) 


-RZ2[1  -f-e'^^R^l-f-Z^  bdawlj-t- 


[l-e-4RL)  2(L2-R2z^) 

2R..  E(L-R) - ^|rg - . 


A(27) 


Carrying  out  the  obvious  cancellation  in  A(27),  followed  by  some  minor  rearrangements,  gives  A(21) 
or  (29). 

DERIVATION  OF  (30): 

Equation  (30)  follows  directly  from  (23),  where  L  =  0,  by  noting  that  aerf(L,  R)  =  2erfR,  by 
multiplying  and  dividing  by  R  after  the  minus  sign,  and  by  using  dxdaw  =  daw(x)/x. 

DERIVATION  OF  (31): 

For  iiic  derivation  of  (31)  use  A(10)  and  A(ll)  to  obtain 

p_  1  f  A(28) 

^  JL-R  '■ 


^  =  e-(w/u)^[2  Lz-z^]^ 


and  make  the  substitution  z  =  i^^  (i  =  ^  -  1 )  to  obtain 
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Using  the  relation 


g2xt-t^  _  y^Hn(x)t° 
^  n!  ’ 


0 


[1,  P.  784] 


and  A(29),  A(28)  becomes 

P  =  — 

With  obvious  modifications  to  A(30)  the  result  for  (31)  follows,  namely, 


P  =  1  fl-e-^'"/ 


_g -(w/u)^[R^-L^]  ^  Hn(  -  I  Lw/u)  2 

VIF  [ 


E 

1 


ni 


where 


Equation  A(31)  is  evaluated  with  the  use  of  recurrence  relations.  Let 
/■  L  -h  R  . 

_  (L-R)"-»E(L-R)-(L-^R)°-^E(L  +  R)  ,  1  ,, 

Vx“  n!  '^2n^*'~2 

=  AT„_2  +  ^5^^^[(L-R)»-l-e-4LR(L  +  R)n-l],  n  >  2, 

To  =  aerf(L,  R),  Tj  =  E(L  -  R)  [l  -  e  '  ^  I'l. 

Let 

Un  =  (^)"(-i)"Hn(-i^L). 

Then  using  the  recurrence  relation  for  Hermite  polynomials  given  in  A(8),  we  have 
Un=[(-i)^(-i)^L]Uo_i  +  2(n  -  l)(f  )2Un  _  2  ,  n  >  2, 

Un  =  2(^)'[-LUn_,+(n-l)Un_2],  n  >  2, 

Uo=l,  Ui=-2(^)2l. 

Therefore  A(31)  can  also  be  written  as 

p=i{(i-.-<»/'')"iR'-‘-'i)«,f(L,R)  Et„ u„l. 


or 


with 


A(29) 


A(30) 


A(31) 


A(32) 


A(33) 
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FORTRAN  LISTINGS  FOR  ELLCOV  AND  SUPPORTING  ROUTINES 


Below  we  give  a  summary  of  the  various  subprograms  that  are  used  for  computing  P.  These 
subprograms  are  listed  in  this  appendix,  except  for  those  which  are  contained  in  NSWCLIB  [13].  The 
NSWCLIB  routines  are  identified  below  with  a  superscript  ♦.  All  subprograms  are  given  in  double 
precision. 


REFERENCING  OF  ROUTINES  USED  TO  COMPUTE  P 


ELLCOV 

uses: 

AERF* 

PKILL* 

CH 

RQUAD 

CIRCV* 

SEQHZ3 

DEPSLN* 

SQUAD 

DPMPAR*  EQSIG 
TQUAl 

EQSIG 

uses: 

AERF* 

DXPARG* 

ERF5 

HSEXP 

FN2 

uses: 

CIRCV* 

PKILL* 

RQUAD 

uses: 

FN2 

SEQHZ3 

uses: 

AERF* 

ERF* 

BDAWl 

ERFCO* 

DAW* 

ERFCR 

DXDAW 

HSEXP* 

DXPARG*  EQSIG 

SQUAD 

uses: 

FN2 

TQUAl 

uses: 

FN2 
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DOUBLE  PRECISION  FUNCTION  ELLCOV(R,HX,HY,HZ,SX,SY,SZ, 

♦  XK0,XK2,SS) 

C- - 

C  (X,Y,Z)  IS  AN  ELEMENT  OF  A  CARTESIAN  COORDINATE  SYSTEM. 

C  ELLCOV  RETURNS  THE  PROBABILITY  OF  A  SHOT  FALLING,  UNDER  AN 

C  ELLIPSOIDAL  NORMAL  DISTRIBUTION,  IN  A  SPHERE  OF  RADIUS  R  WITH 

C  CENTER  (HX,HY,HZ)  AND  RADIUS  R.  THE  DISTRIBUTION  HAS  MEAN 
C  (0,0,0)  AND  STANDARD  DEVIATIONS  SX,  SY  AND  SZ  IN  THE  X,Y,Z 
C  DIRECTIONS,  RESPECTIVELY.  THE  INPUT  PARAMETERS  ARE  R,HX,HY,HZ, 

C  SX,SY,SZ,XK0,XK2,SS,  WHERE  XKO  =  HX*HX-»-HY*HY-l-HZ*HZ, 

C  XK2  =  SQRT(XKO)  AND  SS  =  MAX(SX,SY,SZ).  A  24TH  ORDER  GAUSSIAN 
C  NUMERICAL  INTEGRATION  IS  CARRIED  OUT  IN  RQUAD  AND  SQUAD. 

C 

C  THE  OUTPUT  ELLCOV  IS  ACCURATE  TO  AT  LEAST  6-SIGNIFICANT  DIGITS 
C  WHEN  lD-20  .LT.  ELLCOV  .LT.  .9999999,  AND  (H/SX)**2+(K/SY)**2-|- 
C  (L/SZ)**2  .LT.  2/DPMPAR(l). 

C 

C  REF:  INTEGRATION  OF  THE  TRIVARIATE  NORMAL  DISTRIBUTION  OVER  AN 
C  OFFSET  SPHERE  AND  AN  INVERSE  PROBLEM.  NSWC  TR.  87-27,  2/1988. 

C  REF;  SIGNIFICANT  DIGIT  COMPUTATION  OF  THE  ELLIPSOIDAL  COVERAGE 
C  FUNCTION  AND  ITS  INVERSE.  NAVSWC  TR.  91-487,  8/91. 


EXTERNAL 

EXTERNAL 

« 

AERF  ,CH  ,DEPSLN  ,DPMPAR  ,EQSIG  ,PKILL 
RQUAD  ,SEQHZ3  ,SQUAD  ,TQUA1  ,CIRCV 

DOUBLE  PRECISION 

AERF 

,DEPSLN 

,DPMPAR  ,EQSIG 

DOUBLE  PRECISION 

X(12) 

,Y(12) 

DOUBLE  PRECISION 

K 

,L 

DOUBLE  PRECISION 

A 

,AH 

,AK 

,AL 

,A2 

,A3 

DOUBLE  PRECISION 

B 

,BH 

,BK 

,BL 

,B1 

,C5 

DOUBLE  PRECISION 

C6 

,C7 

,C8 

,C9 

,D 

,DEP 

DOUBLE  PRECISION 

EPS 

,EI 

,FH 

,FK 

,FL 

,F2 

DOUBLE  PRECISION 

H 

,HH 

,HK 

,HL 

,HX 

,HY 

DOUBLE  PRECISION 

HZ 

,PY 

,R 

,R1 

,R2 

,R3 

DOUBLE  PRECISION 

SS 

,SX 

,SY 

,SZ 

,S1 

,S2 

DOUBLE  PRECISION 

S3 

,T5 

,T8 

,W1 

,XK0 

,XK2 

DOUBLE  PRECISION 

XK2 

,XK7 

,ZH 

,ZL 

,Z4 

X(*),Y{*)  -  GAUSSIAN  ABSCISSAS  AND  WEIGHTS  OF  ORDER  24  ON  (-1,1). 


♦ 

♦ 

* 

* 

* 


C^ 


DATA  X(l)  /6.40568928626056D-2/  ,X(2)  /.191118867473616D0/, 
X(3)  /.31 5042679696 163D0/  ,X(4)  /.433793507626045D0/, 
X(5)  /.545421471388840D0/  ,X(6)  /.648093651936976DO/, 
X(7)  /.740124191578554D0/  ,X(8)  /.820001985973903D0/, 
X(9)  /.886415527004401D0/  ,X(10)  /.938274552002733D0/, 
X(ll)  /.974728555971309D0/  ,X(12)  /.995187219997021D0/ 
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DATA  Y(l)  /.127938195346752D0/  ,Y(2)  /.125837456346828D0/, 

*  Y(3)  /.121670472927803D0/  ,Y(4)  /.115505668053726D0/, 

*  Y(5)  /.  1074442701 15966D0/  ,Y(6)  /9.76186521041139D-2/, 

*  Y(7)  /8.61901615319533D-2/  ,Y(8)  /7.33464814110803D-2/, 

*  Y(9)  /5.92985849154368D-2/  ,Y(10)  /4.42774388174198D-2/, 

*  Y(ll)  /2.85313886289337D-2/  ,Y(12)  /1.23412297999872D-2/ 


DATA  B1  /0.5641895835477563D0/,C5  /7.0D0/, 

♦  C6  /0.3989422804014327D0/,C7  /0.7071O67811865475D0/, 

•  C8  /1.414213562373095D0/  ,C9  /9.6D0/,  F2  /IDl/ 


B1  =  1/SQRT(PI),  C6  =  1/SQRT(2PI) 


EPS  =  10*DPMPAR(1) 

T8  =  1D2*DPMPAR(2) 

Z4  =  MAX(T8,lD-25*lD-25) 
DEP  =  -DEPSLN(O) 

N  =  12 
J1  =  0 

H  =  ABS(HX) 

K  =  ABS(HY) 

L  =  ABS(HZ) 

51  =  SX 

52  =  SY 

53  =  SZ 
PY  =  O.DO 
ELLCOV  =  O.ODO 

El  =  1.5D0*SX*SY*SZ*Z4/C6 
XK7  =  R*R*R 
IF  (XK7  .LE.  El)  RETURN 
A3  =  (XK2  -  R)*C7/SS 
A  =  EXP(-A3*A3) 

IF  (XK2  .LT.  R)  GO  TO  5 
IF  (XK7*A  .LE.  El)  RETURN 
GO  TO  10 


5  PY  =  EQSIG(R,XK0,XK2,SS,EPS,DEP,A) 
El  =  DMAX1(1D-11,  5*EPS) 

IF  (PY  .LT.  l.ODO  -  El)  GO  TO  10 

ELLCOV  =  l.ODO 

RETURN 

10  IF  (SX  .NE.  SY  .OR.  SY  .NE.  SZ)  GOTO  15 


SI  =  S2  =  S3. 


ELLCOV  =  PY 

IF  (R  .GT.  XK2)  RETURN 

ELLCOV  =  EQSIG(R,XK0,XK2,SS,EPS,DEP,A) 

RETURN 
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C - 

C  SMALL  R 

C - 

15  R1  =  R/Sl 
R2  =  R/S2 
R3  =  R/S3 
HR  =  H/Sl 
HK  =  K/S2 
HL  =  L/S3 
AH  =  RURl 
BH  =  HH*HH 
AK  =  R2*R2 
BK  =  HK*HK 
AL  =  R3*R3 
BL  =  HL*HL 

T5  =  AH*(BH  -  IDO)  +  AK*(BK  -  IDO)  +  AL*(BL  -  IDO) 

IF  (ABS(T5)  .GT.  lD-3)  GOTO  20 

W1  =  4*(AH*AH*(BH  -  .5D0)  -|-  AK*AK*(BK  -  .5D0)  + 

♦  AL*AL*(BL  -  .5D0)) 

AH  =  (T5*T5  +  ABS(Wl))/280 

IF  (AH  .GT.  MAX(EPS,2.5D-8))  GOTO  20 

ELLCOV  =  B1»R1*R2*R3/(3*C7)*(1.D0  -1-  T5/10  +  (T5*T5  -  Wl)/280)* 

♦  EXP(-(BH  -I-  BK  -t-  BL)/2) 

RETURN 

20  Jl  =  3 

IF  (SX  .EQ.  SY  .AND.  H  4-  K  .EQ.  O.DO)  GOTO  25 
Jl  =  2 

IF  (SX  .EQ.  SZ  .AND.  H  -t-  L  .EQ.  O.DO)  GOTO  25 
Jl  =  1 

IF  (SZ  .NE.  SY  .OR.  L  4-  K  .NE.  O.DO)  GOTO  30 
25  CALL  CH(J1,H,K,L,S1,S2,S3) 

C- - 

C  SI  =  S2  AND  H  =  K  =  0. 

C - 

CALL  SEQHZ3(R,L,S1,S3, EPS, DEP,XK0,A,PY, ELLCOV, KKK) 

RETURN 

C- - 

30  J  =  2 
II  =  0 

AH  =  MAX(-C5,(HH  -  R1)*C7) 

IF  (AH  .GT.  C9)  RETURN 
BH  =  MIN(F2,HH*C7) 

AK  =  MAX(-C5,(HK  -  R2)*C7) 

IF  (AK  .GT.  C9)  RETURN 
BK  =  MIN(F2,HK*C7) 

AL  =  MAX(-C5,(HL  -  R3)*C7) 

IF  (AL  .GT.  C9)  RETURN 
BL  =  MIN(F2,HL*C7) 


B-5 


NAVSWC  TR  91-487 


T8  =  MAX(1D-42,T8) 

ZL  =  MIN(EPS/30,1D-12) 

CALL  TQUA1(AH,BH,N,X,R,K,L,H,S2,S3,S1,PY,T8,DEP,FH) 
IF  (AH  .GT.  -6D0)  GOTO  35 
IF  ((HH  -  C8*AH)+*2  .GT.  R1*R1*ZL)  GOTO  35 
II  =  1 
GOTO  110 
35  T5  =  BH  -  AH 
N1  =  0 

IF  (T5  .NE.  O.DO)  GOTO  40 
N1  =  1 
GOTO  45 
40  J1  =  1 
A  =  AH 
B  =  BH 

45  CALL  TQUA1(AK,BK,N,X,R,L,H,K,S3,S1,S2,PY,T8,DEP,FK) 
IF  (AK  .GT.  -6D0)  GOTO  50 
IF  ((HK  -  C8+AK)**2  .GT.  R2*R2*ZL)  GOTO  50 
II  =  1 
J1  =  2 
GOTO  no 
50  A2  =  BK  -  AK 
N2  =  0 

IF  (A2  .NF  u.DO)  GOTO  55 
N2  =  1 
GOTO  85 

55  IF  (N1  .NE.  0)  GOTO  65 

If  (AH  .LT.  -2D0  .AND.  AK  -  AH  .GT.  4D0)  GOTO  85 
IF  (AK  .LT.  -2D0  .AND.  AH  -  AK  .GT.  4D0)  GOTO  65 
W1  =  T5  -  A2 

IF  (ABS(Wl)  .GT.  5D2*EPS)  GOTO  60 
IF  (AH  .GT.  AK  .AND.  AH  .LT.  -2.D0)  GOTO  85 
60  IF  (Wl)  65,65,85 
65  J1  =  2 
A  =  AK 
B  =  BK 

CALL  TQUA1(AL,BL,N,X,R,H,K,L,S1,S2,S3,PY,T8,DEP,FL) 

IF  (AL  .GT.  -6D0)  GOTO  70 

IF  ((HL  -  C8*AL)**2  .GT.  R3*R3*ZL)  GOTO  70 

II  =  1 

J1  =  3 

GOTO  no 

70  T5  =  BL  -  AL 

IF  (T5  .NE.  O.DO)  GOTO  75 
GOTO  no 

75  IF  (N2  .NE.  0)  GOTO  105 

IF  (AK  .LT.  -2D0  .AND.  AL  -  AK  .GT.  4D0)  GOTO  110 
IF  (AL  .LT.  -2D0  .AND.  AK  -  AL  .GT.  4D0)  GOTO  105 
Wl  =  T5  -  A2 
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IF  (ABS(Wl)  .GT.  5D2*EPS)  GOTO  80 
IF  (AK  .GT.  AL  .AND.  AK  .LT.  -2.D0)  GOTO  110 
80  IF  (Wl)  110,105,105 

85  CALL  TQUA1(AL,BL,N,X,R,H,K,L,S1,S2,S3,PY,T8,DEP,FL) 
IF  (AL  .GT.  -6D0)  GOTO  90 
IF  ((HL  -  C8*AL)**2  .GT.  R3*R3*ZL)  GOTO  90 
II  =  1 
J1  =  3 
GOTO  110 
90  A2  =  BL  -  AL 

IF  (A2  .NE.  O.DO)  GOTO  95 
IF  (N1  .NE.  1)  GOTO  110 
ELLCOV  =  O.DO 
RETURN 

95  IF  (N1  .EQ.  1)  GOTO  105 

IF  (AH  .LT.  -2D0  .AND.  AL  -  AH  .GT.  4D0)  GOTO  110 
IF  (AL  .LT.  -2D0  .AND.  AH  -  AL  .GT.  4D0)  GOTO  105 
Wl  =  T5  -  A2 

IF  (ABS(Wl)  .GT.  5D2*EPS)  GOTO  100 
IF  (AH  .GT.  AL  .AND.  AH  .LT.  -2.D0)  GOTO  110 
100  IF  (Wl)  105,105,110 
105  J1  =  3 
A  =  AL 
B  =  BL 

C- - 

no  CALL  CH(J1,H,K,L,S1,S2,S3) 

IF  (SI  .NE.  S2)  GOTO  115 

C - 

C  SI  =  S2,  J  =  1. 


J  =  1 

D  =  SQRT(H*H  -I-  K*K)/S1 
GO  TO  125 

115  IF  (H -1- K  .NE.  O.DO)  GOTO  125 

C  H  =  K  =  0,  J  =  0. 

IF  (SI  .GT.  S2)  GO  TO  120 
Wl  =  SI 

51  =  S2 

52  =  Wl 
120  D  =  S2/S1 

J  =  0 

125  IF  (II  .EQ.  0)  GOTO  135 

C  SEE  P.IO  OF  REPT.  87-27. 

C - 

IF  (J  .LT.  2)  GO  TO  130 

CALL  PKILL(R,Sl,S2,H,K,ELLCOV) 
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RETURN 

130  CALL  C1RCV(R/S1,D,J,ELLC0V,IER) 

RETURN 

C - 

C  PKILL  OR  CIRCV  .GT.  1  -  lE-8  FOR  ALL  T  .GE.  T8. 

C - 

135  PY  =  O.ODO 

XK7  =  DMAX1(S1,S2) 

T8  =  C5*XK7  -I-  SQRT(H*H  +  K*K) 

T8  =  R*R  -  T8*T8 
ZL  =  L/(C8+S3) 

IF  (T8  .LT.  O.ODO)  GO  TO  140 
A3  =  SQRT(T8) 

ZH  =  A3/(C8*S3) 

T8  =  ZL  -  ZH 

IF  (T8  .GE.  B)  GO  TO  140 

B  =  T8 

PY  =  AERF(ZL,ZH)/2 
ELLCOV  =  PY 

IF  ((l.DO  -  PY)  .LT.  El)  RETURN 

IF  (B  -  A  .LE.  ABS(B  -(-  A)*DMAX1(1D-9,EPS))  RETURN 

C  OBTAIN  SHARP  INTEGRATION  LIMITS  FOR  A  AND  B. 

C- - 

140  A2  =  B  -  A 
XK7  =  A 
T8  =  B 

CALL  SQUAD(A,B,W1,N,X,Y,R,H,K,L,S1,S2,S3,EPS,DEP, 

♦  J,D,PY,1L) 

IF  (W1  .LT.  Z4)  GOTO  145 

C  IF  A  AND  B  ARE  UNCHANGED  BY  SQUAD,  ELLCOV  =  W1  +  PY. 

IF(ABS(B  -  A)  .LT.  A2*(1.D0  -  EPS))  GOTO  150 
145  ELLCOV  =  PY  -(-  W1 
RETURN 

C  START  OF  GAUSSIAN  INTEGRATION. 

C - 

C  II  .GE.  1.  GIVES  NO.  OF  SUBINTERVALS  OVER  WHICH  GAUSSIAN 
C  NUMERICAL  INTEGRATION  IS  APPLIED.  GENERALLY  II  =  1.  HIGHER 
C  VALUES  USED  ONLY  FOR  CHECKING  PURPOSES. 

C - 

150  II  =  1 

CALL  RQUAD(A,B,ELLC0V,I1,N,X,Y,R,H,K,L,S1,S2,S3,DEP,J,D) 

ELLCOV  =  BUELLCOV  +  PY 

IF  (ELLCOV  .GT.  l.ODO)  ELLCOV  =  l.ODO 

RETURN 

END 
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DOUBLE  PRECISION  FUNCTION  BDAW1(XM,XP,DEL,T) 

C  LET  FI  =  XP/DEL,  SO  =  XM/DEL. 

C  DAW(X)  =  EXP(-X*X)  ♦  INTEGRAL  (FROM  0  TO  X)  EXP(W*V')  DW.  BDAWl 
C  GIVES  THE  QUANTITY  (XM/XP**2)  FRAC(XP)  -  T*(XP/XM**2)  FRAC(XM), 

C  WHERE  FRAC(X)  APPEARS  IN  THE  MINIMAX  RATIONAL  APPROXIMATION 
C  GIVEN  FOR  DAW(X)  BY  (1/(2X)  -1-  1/(2  X**3))  FRAC(X),  WHERE  X  .GE.  5. 

C  BDAWl  IS  CALLED  FROM  SEQHZ3,  WHICH  IS  CALLED  BY  ELLCOV. 

C  EVALUATION  OF  DAWSON’S  INTEGRAL  FOR  ALL  REAL  X  BY  DAW(X)IS  BASED 
C  ON  RATIONAL  CHEBYSHEV  APPROXIMATIONS  PUBLISHED  IN  MATH.  COMP. 

C  24,  171-178(1970)  BY  CODY,  PACIOREK  AND  THACHER. 


DIMENSION  P4(7)  ,Q4(6) 

C- - 

DOUBLE  PRECISION  DEL  ,FRM  ,FRP  ,P4  ,Q4  ,T 

DOUBLE  PRECISION  XLG  ,XM  ,XP  ,ZM2  ,ZP2 


DATA  XLG/16777216.0D0/ 


COEFFICIENTS  FOR  R(6,6)  APPROXIMATION, 

IN  J-FRACTION  FORM,  USED  FOR  ABS(X)  .GT.  5.0 


DATA  P4(l)/-.315576735766984D-l-02/,  P4(2)/-.10079I496592972D+02/, 

*  P4(3)/-.7107I3709224200D-(-01/,  P4(4)/-.596879853243925D4-01/, 

*  P4(5)/-.449773645376092D-1-01/,  P4(6)/-.249999965398199D4-01/, 

*  P4(7)/  .499999999999330D-f00/ 

DATA  Q4(l)/  .168874162I55616D-F03/,  Q4(2)/  .698280748271071D-(-01/, 

*  Q4(3)/-.2I302962I139I81D4-02/,  Q4(4)/-.712I57348463305D-t-01/, 

*  Q4(5)/-.250005973I92356D-f-01/,  Q4(6)/  .750000000715687D-I-00/ 


FRP  =  O.ODO 
BDAWl  =  O.ODO 

IF  (DEL*XLG  .LE.  ABS(XM))  RETURN 
IF  (DEL+XLG  .LE.  XP)  GOTO  10 
ZP2  =  (DEL/XP)**2 
DO  5  I  =  1,  6 

5  FRP  =  ZP2*Q4(I)  /  (ZP2*(P4(I)  -|-  FRP)  -f-  1) 
FRP  =  P4(7)  -I-  FRP 

10  FRM  =  O.DO 
IF  (T  .EQ.  O.DO)  GOTO  25 
ZM2  =  (DEL/XM)**2 
DO  15  I  =  I,  6 

15  FRM  =  ZM2+Q4(I)  /  (ZM2*(P4(I)  -I-  FRM)  -h  1) 
FRM  =  P4(7)  -I-  FRM 

20  BDAWl  =  XM*FRP/XP**2  -  T*XP*FRM/XM**2 
RETURN 

25  BDAWl  =  XM*FRP/XP**2 
RETURN 
END 
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SUBROUTINE  CH(J,H,K,L,S1,S2,S3) 


BASED  ON  THE  VALUE  OF  J  (=1,2,3)  H,K,L  AND  S1,S2,S3  ARE 
INTERCHANGED.  USED  IN  ELLCOV. 


DOUBLE  PRECISION  K  ,L 

C - 

DOUBLE  PRECISION  H  ,S1  ,S2  ,S3  ,X 

O - 

IF  (J  .EQ.  3)  RETURN 
IF  (J  .EQ.  2)  GO  TO  5 
X  =  H 
H  =  L 
L  =  X 
X  =  SI 

51  =  S3 
S3  =  X 
RETURN 

5  X  =  K 
K  =  L 
L  =  X 
X  =  S2 

52  =  S3 

53  =  X 

10  RETURN 
END 
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DOUBLE  PRECISION  FUNCTION  DXDAW(X) 

C - 

C  THIS  FUNCTION  COMPUTES  VALUES  OF  THE  FUNCTION  - 
C  EXP(-X*X)  *  (1/X)*INTEGRAL  (FROM  0  TO  X)  EXP(T*T)  DT, 

C  DEFINED  FOR  ALL  REAL  ARGUMENTS.  USED  IN  SEQHZ3. 

C  THE  MAIN  COMPUTATION  INVOLVES  EVALUATION  OF  RATIONAL  CHEBYSHEV 
C  APPROXIMATIONS  PUBLISHED  IN  MATH.  COMP.  24,  171-178(1970)  BY 
C  CODY,  PACIOREK  AND  THACHER. 

C- - 

DOUBLE  PRECISION  Pl(9)  ,Q1(9)  ,P2(8)  ,Q2(7)  ,P3(8)  ,Q3(7)  ,P4(7)  ,Q4(6) 

DOUBLE  PRECISION  FRAC  ,SUMP  ,SUMQ  ,W2  ,X  ,Y  ,XLARGE 

C - 

DATA  XLARGE/16777216.0D0/ 


COEFFICIENTS  FOR  R(8,8)  APPROXIMATION, 
USED  FOR  ABS(X)  .LT.  2.5 


DATA 


DATA 


♦ 


P 1  ( 1 )/.  1  OOOOOOOOOOOOOOD-l-01/, 

P 1  (3)/.456738974064825D-0 1  /, 

Pl(5)/.360079463580992D-03/, 

P 1 (7)/.634674256878843D-06/, 

P1(9)/.977985913592343D-10/ 

Q1(1)/.100000000000000D+01/, 

Q  1(3)/.  133052308640737D-f  00/, 

Ql(5)/.220437428972266D-02/, 

Ql(7)/.887964712053131D-05/, 

Ql(9)/.574807177698046D-08/ 


Pl(2)/-.135599049815353D-t-00/, 

Pl(4)/-.258323495918051D-02/, 

Pl(6)/-.944375029163387D-05/, 

Pl(8)/-.711645839183817D-08/, 

Q 1  (2)/. 53 1 0676 168513 1 0D-(-00/, 
Q1(4)/.206907491644210D-01/, 
Q1  (6)/.  16670680 1 664365D-03/, 
Ql(8)/.311750854173480D-06/, 


COEFFICIENTS  FOR  R(7,7)  APPROX WATION, 
IN  J-FRACTION  FORM,  USED  FOR 
2.5  .LE.  ABS(X)  .LT.  3.5 


DATA  P2(l)/-.150695651187161D-(-01/,  P2(2)/  .293365747395449D-I-02/, 

♦  P2(3)/-.400000893643550D+02/,  P2(4)/-.757931918089369D-01/, 

♦  P2(5)/-.889106479747812D-»-01/,  P2(6)/  .152644099623699D-t-02A 

♦  P2(7)/-.597678086823489D-f01/,  P2(8)/  .500236896088668D-1-00/ 

DATA  Q2(l)/-.673106069744813D-»-00/,  Q2(2)/  .124486788262252D-I-04/, 

♦  Q2(3)/  .721193217600229D-I-01/,  Q2(4)/  .1124616620245750^-03/, 

♦  Q2(5)/  .7291 775564 15532D-t-02/,  Q2(6)/  .115840292551888D-f-03/, 

♦  Q2(7)/  .226064666074309D+00/ 


COEFFICIENTS  FOR  R(7,7)  APPROXIMATION, 
IN  J-FRACTION  FORM,  USED  FOR 
3.5  .LE.  ABS(X)  .LE.  5.0 
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DATA  P3(l)/  .4764056452732290-1-01/,  P3(2)/-.266167674896399D-^02/, 

♦  P3(3)/-.916804879813552D-t-01/,  P3(4)/-.150507703496692D-l-02/, 

♦  P3(5)/  .506460153742231D+01/,  P3(6)/-.498544802986608D-f01/, 

♦  P3(7)/-.149838042036691D-t-01/,  P3(8)/  .4999999027050540-1-00/ 

DATA  Q3(l)/  .2877761229731870+03/,  Q3(2)/  .2561057223422260+02/, 

♦  Q3(3)/  .7517012777440670+02/,  Q3(4)/  .1465151677831090+03/, 

♦  Q3(5)/  .3307077246761140+02/,  Q3(6)/-.14871581 17871950+01/, 

♦  Q3(7)/  .2500114596118390+00/ 


COEFFICIENTS  FOR  R(6,6)  APPROXIMATION, 

IN  J-FRACTION  FORM,  USED  FOR  ABS(X)  .GT.  5.0 


DATA  P4(  1  )/-.3 155767357669840+02/,  P4(2)/-.  1007914965929720+02/, 

♦  P4(3)/-.710713709224200D+01/,  P4(4)/-.596879853243925D+01/, 

♦  P4(5)/-.449773645376092D+01/,  P4(6)/-.249999965398199D+01/, 

♦  P4(7)/  .4999999999993300+00/ 

DATA  Q4(l)/  .1688741621556160+03/,  Q4(2)/  .6982807482710710+01/, 

♦  y4(3)/-.213029621139181D+02/,  Q4(4)/-.712157348463305D+01/, 

♦  Q4(5)/-.250005973192356D+01/,  Q4(6)/  .7500000007156870+00/ 


Y  =  X  *  X 

IF  (ABS(X)  .GT.  XLARGE)  GO  TO  35 
IF  (Y  .GE.  6.2500)  GO  TO  5 


ABS(X)  .LT.  2.5 

SUMP  =  ((({(((Pl(9)  ♦  Y  +  Pl(8))  *  Y  +  Pl(7))  *  Y  +  Pl(6)) 

1  ♦  Y  +  Pl(5))  *  Y  +  Pl(4))  *  Y  +  Pl(3))  ♦  Y  +  Pl(2)) 

2  ♦Y  +  Plil) 

SUMQ  =  (((((((Ql(9)  ♦  Y  +  Ql(8))  ♦  Y  +  Ql(7))  *  Y  +  Ql(6)) 

1  *  Y  +  Ql(5))  *  Y  +  Ql(4))  ♦  Y  +  Ql(3))  *  Y  +  Ql(2)) 

2  ♦  Y  +  Ql(l) 

DXDAW  =  SUMP  /  SUMQ 
RETURN 


2.5  .LE.  ABS(X)  .LT.  3.5 


5  IF  (Y  .GE.  12.2500)  GO  TO  15 
FRAC  =  0.000 
DO  10  I  =  1,  7 

10  FRAC  =  Q2(I)  /  (P2(I)  +  Y  +  FRAC) 
DXDAW  =  (P2(8)  +  FRAC)  /  Y 
RETURN 


3.5  .LE.  ABS(X)  .LT.  5.0 
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15  IF  (Y  .GE.  25.0D0)  GO  TO  25 
FRAG  =  O.ODO 
DO  20  I  =  1,  7 

20  FRAG  =  Q3(I)  /  (P3(I)  -I-  Y  -I-  FRAG) 
DXDAW  =  (P3(8)  +  FRAG)  /  Y 
RETURN 


5.0  .LE.  ABS(X)  .LE.  XLARGE 


25  W2  =  l.ODO  /  X  /  X 
FRAG  =  O.ODO 
DO  30  I  =  1,  6 

30  FRAG  =  Q4(I)  /  (P4(I)  +  \  +  FRAG) 

FRAG  =  P4(7)  -t-  FRAG 

DXDAW  =  (0.5D0  +  0.5D0  *  W2  ♦  FRAG)  /  Y 

RETURN 


XLARGE  .LT.  ABS(X) 


35  DXDAW  =  0.5D0  /  Y 
RETURN 
END 
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DOUBLE  PRECISION  FUNCTION  EQSIG(R3,XK0,XK2,S,EPS,DEP,A5) 

C - 

C  EQSIG  COMPUTES  ELLCOV  FOR  EQUAL  SIGMAS.  XKO  =  H*H-(-K*K-|-L+L, 
C  XK2  =  SQRT(XKO).  EPS  =  10*DPMPAR(1). 

C  A  =  (XK2-R3)*SQRT(.5)/S,  A5  =  EXP(-A*A). 

C - 

EXTERNAL  ERF5  ,HSEXP  ,AERF  .DXPARG 
C - 


DOUBLE  PRECISION  AERF  , DXPARG 


DOUBLE  PRECISION 

A3 

,A4 

,A5 

,A6 

,C7 

,DEP 

DOUBLE  PRECISION 

E 

,EPS 

,G1 

,H2 

,H4 

,H5 

DOUBLE  PRECISION 

H6 

,H7 

,P4 

,P5 

,Q5 

,R3 

DOUBLE  PRECISION 

S 

,si 

,T 

,V1 

,XK0 

,XK2 

DOUBLE  PRECISION 

z 

DATA  A4/1.1283  79167  09551  257D0/,  C7/.70710  67811  86547  524D0/ 


A4  =  2/SQRT(PI) 


EQSIG  =  O.DO 
H6  =  R3*C7/S 
A6  =  XK2*C7/S 
A3  =  A6  -  H6 
E  =  MAX(1D-11,EPS) 

IF  (XKO  .NE.  O.ODO)  GOTO  15 
IF  (H6  .GT.  .071D0)  GOTO  10 


D  =  0  AND  R3/(SQRT(2)*S)  .LE.  .071 


N  =  3 

EQSIG  =  l.DO 
A6  =  l.DO 
P5  =  2.D0*H6*H6 
Q5  =  P5*H6/N 
5  N  =  N  4-  2 

A6  =  A6*P5/N 
EQSIG  =  EQSIG  +  A6 
IF  (A6  .GT.  E+EQSIG)  GOTO  5 
EQSIG  =  A4*A5i.Q5*EQSIG 
KKK  =  13 
RETURN 

D  =  0  AND  R3/(SQRT(2)*S)  .GT.  .071 


10  CALL  ERF5(-A3,A5,P4) 
EQSIG  =  P4  -  A4*H6*A5 
KKK  =  14 
RETURN 
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D  .NE.  0  AND  R3/(SQRT(2)*S)  .GE.  .425 


15  IF  (H6  .LT.  .425D0)  GOTO  30 
P5  =  4.0D0*A6*H6 
IF  (A3  .LT.  4.5D0)  GOTO  25 
IF  (P5  .LT.  DEP)  GOTO  25 
SI  =  H6/A6 
E  =  MAX(2.5D-8,EPS) 

Z  =  l.DO 
N  =  -1 

H7  =  2*A3*A3 
20  N  =  N  4-  2 
Z  = -N*Z/H7 
SI  =  SI  -t-  Z 

IF  (ABS(Z)  .GT.  E*ABS(S1))  GOTO  20 
SI  =  SI  -  .5D0*Z 
EQSIG  =  .25D0*A4*A5/A3*S1 
RETURN 

25  P4  =  AERF(A6,H6) 

VI  =  .1.D0 

CALL  HSEXP(-P5,V1,Q5) 

EQSIG  =  0.5DO*P4  -  A4»H6*A5*Q5 

KKK  =  15 

RETURN 


D  .NE.  0  AND  R3/(SQRT(2)*S)  .LT.  .425 


30  T  =  A6*A6/2 
KKK  =  16 

IF  (T  4-  10  .GT.  -DXPARG(l))  RETURN 
G1  =  EXP(-T) 

H2  =  H6*H6 
H7  =  H2*Hfi 
H5  =  2*G1 
H4  =  G1 
N  =  1 

VI  =  1.5D0 
SI  =  H5*H7/V1 
T  =  2*T 
J  =  0 

35  H7  =  H2*H7 

H4  =  (T*H5  -  H4)/N 
H5  =  (H4  -  H5)/V1 
N  =  N  4-  1 
VI  =  N  4-  .5D0 
Z  =  H5*H7/V1 
SI  =  SI  4-  Z 
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IF  (ABS(Z)  .GT.  E+ABS(S1))  GOTO  35 
IF  (J  .GT.  0)  GOTO  40 
J  =  1 
GOTO  35 

40  EQSIG  =  A4*Gl*Sl/2 
KKK  =  17 
RETURN 
END 
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SUBROUTINE  ERF5(X,EXPP,Y) 


C  Y  =  THE  REAL  ERROR  FUNCTION  OF  X.  EXPP  =  EXP(-X*X). 
C  IF  EXP(-X*X)  NOT  AVAILABLE  SET  EXPP  .LT.  0. 


P— . 

DIMENSION 

A(4) 

.B(4) 

,P(8) 

.Q(8) 

,R(5) 

,S(5) 

DOUBLE  PRECISION 

A 

,AX 

.B 

,BOT 

,c 

,EXPP 

DOUBLE  PRECISION 

P 

,Q 

,R 

,s 

,T 

,TOP 

c~ 

DOUBLE  PRECISION 

X 

,X2 

,Y 

DATA  C/.564189583547756D0/ 

DATA  A(l)/-1.65581836870402D-4/,  A(2)/3.25324098357738D-2/, 

*  A(3)/1.02201136918406D-1/.  A(4)/1.12837916709552D0/ 

DATA  B(l)/4.64988945913179D-3/,  B(2)/7.01333417158511D-2/, 

*  B(3)/4.23906732683201D-1/,  B(4)/1.00000000000000D0/ 

DATA  P(l)/-1.36864857382717D-7/,  P(2)/5.64195517478974D-1/, 

*  P(3)/7.21175825088309D0/,  P(4)/4.31622272220567D01/, 

»  P(5)/1.52989285046940D02/,  P(6)/3.39320816734344D02/, 

*  P(7)/4.51918953711873D02/,  P(8)/3.00459261020162D02/ 

DATA  Q(1)/1.00000000000000DO/,  Q(2)/1.27827273196294D01/, 

*  Q{3)/7.70001529352295D01/,  Q(4)/2.77585444743988D02/, 

*  Q(5)/6.38980264465631D02/,  Q(6)/9.31354094850610D02/, 

*  Q(7)/7.90950925327898D02/,  Q(8)/3.00459260956983D02/ 

DATA  R(1)/2.10144126479064D0/,  R{2)/2.62370141675169D01/, 

*  R(3)/2.13688200555087D01/,  R(4)/4.65807828718470D0/, 

*  R(5)/2.82094791773523D-1/ 

DATA  S(1)/9.41537750555460D01/,  S(2)/1.87114811799590D02/, 

*  S(3)/9.90191814623914D01/,  S(4)/1.80124575948747D01/, 

*  S(5)/1.00000000000000DO/ 

c - 

AX  =  ABS(X) 

T  =  X*X 

IF  (AX  .GE.  0.5D0)  GO  TO  5 

TOP  =  ({A(1)*T  -(-  A(2))*T  +  A(3))*T  +  A(4) 

BOT  =  ((B(1)*T  +  B(2))*T  -I-  B(3))*T  +  B(4) 

Y  =  X*TOP/BOT 
RETURN 

C 

5  IF  (AX  .GT.  4.0D0)  GO  TO  10 

TOP  =  ((((((P(1)*AX  +  P(2))*AX  +  P(3))*AX  +  P(4))*AX  +  P(5))*AX 

*  +  P(6))*AX  +  P(7))*AX  -I-  P(8) 

BOT  =  ((((((Q(1)*AX  Q(2))*AX  H-  Q(3))*AX  -I-  Q(4))*AX  4-  Q(5))*AX 

*  4-  Q(6))*AX  4-  Q(7))*AX  4-  Q(8) 

IF  (EXPP  .LT.  O.DO)  EXPP  =  EXP(-T) 

Y  =  0.5D0  4-  (0.5D0  -  EXPP*TOP/BOT) 

IF  (X  .LT.  O.ODO)  Y  =  -Y 

RETURN 
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10  Y  =  l.ODO 

IF  (AX  .GE.  5.6D0)  GO  TO  15 
X2  =  l.ODO/T 

TOP  =  (((R(1)*X2  -I-  R(2))*X2  +  R(3))*X2  +  R(4))*X2  +  R(5) 
BOT  =  (((S(1)*X2  -I-  S(2))*X2  +  S(3))*X2  +  S(4))*X2  -I-  S(5) 

Y  =  (C  -  TOP/(T*BOT))  /  AX 

IF  (EXPP  .LT.  O.DO)  EXPP  =  EXP(-T) 

Y  =  0.5D0  -I-  (0.5D0  -  EXPP*Y) 

15  IF  (X  .LT.  O.ODO)  Y  =  -Y 

RETURN 

END 
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DOUBLE  PRECISION  FUNCTION  FN2(T5,1M,IL,R3,H,K,L, 

♦  S1,S2,S3,HZ,DEP,J,D) 

C- - 

C  FN2  GIVES  THE  INTEGRAND  OF  ELLCOV.  USED  BY  RQUAD, SQUAD, TQUADl. 
C  J  =  0,1,2  SPECIFIES  CIRCV  OR  PKILL.  DEP  =  -DEPSLN(O).  D  NEEDED 
C  FOR  CIRCV.  HZ  =  L*SQRT(l/2)/S3.  SEE  BELOW  FOR  IM,  IL. 

C - 

C  U(J)  ARE  USED  WHEN  SQUAD  CALLS  FN2  AND  THE  INPUT  INTEGRATION 
C  LIMITS  A,B  ARE  UNCHANGED  FROM  L-R  AND  L.THEN  R,  USED  IN  PKILL 
C  OR  CIRCV,  IS  GIVEN  BY  R3*U(IM). 

C  IF  IL  .NE.  0  COMPUTE  FN2  USING  THE  U(J)  USING  INPUT  IM. 

C  IF  IL  .EQ.  0  COMPUTE  FN2  USING  T5  WITHOUT  USING  U(J). 

C  U(I)=  SQRT(1  -  .25*(1  +  X(I))**2),  I  =  12 . 1. 

C  U(N-(-I)=  SQRT(1  -  .25*(1  -  X(I))**2), 

C  X(I)  ARE  THE  GAUSSIAN  ABSCISSAS  OF  ORDER  24.  N  =  12. 


EXTERNAL 

DOUBLE  PRECISION 

CIRCV  ,PKILL 
K  ,L 

,U(24) 

DOUBLE  PRECISION 

C8 

,D 

,DEP 

,HZ 

,P 

DOUBLE  PRECISION 

R 

,R3 

,S1 

,S2 

,S3 

,T5 

DOUBLE  PRECISION 

V2 

,XJI 

C - 

DATAU(12)  /.6933245481114434D-1/,  U(ll) /.1584669762375325D0/, 

♦  U(10) /.2465216831531283D0/ ,  U(9)  /.3322034239275342D0/, 

♦  U(8)  /.4146060693752110D0/  ,  U(7)  /.4929421360269080D0/, 

♦  U(6)  /.5665216929750006D0/ ,  U(5)  /.6347583153788445D0/, 

♦  U(4)  /.6971793487850529D0/ ,  U(3)  /.7534359213923792D0/, 

♦  U(2)  /.8033112478280709D0/  ,  U(l)  /.8467264801504051D0/, 

♦  U(13)  /.8837435290006371DO/  ,  U(14)  /.9145642833397272D0/, 

♦  U(15)  /.9395256076023404DO/  ,  U(16)  /.9590894389984678D0/, 

♦  U(17)  /.9738272897857341D0/  ,  U(18)  /.9843985374573839D0/, 

♦  U(19) /.9915221334137353D0/ ,  U(20) /.9959418550983111D0/, 

♦  U(21)  /.9983860184685972D0/  ,  U(22)  /.9995236326707759D0/, 

♦  U(23)  /.9999201660778605D0/  ,  U(24)  /.9999971046393888D0/ 

C- - 

DATA  C8  /1.4142  13562  37310DO/ 

C - 

FN2  =  O.DO 

IF  (T5  .GT.  11.2D0)  RETURN 
XJl  =  EXP(-T5*T5) 

P  =  HZ  -  T5 

IF  (IL  .EQ.  0)  GOTO  5 

R  =  R3*U(IM) 

GOTO  10 
5  V2  =  C8»S3*T5 

R  =  (R3  -  L)*(R3  -h  L)  -f-  V2*(2*L  -  V2) 


B-19 


NAVSWC  TR  91-487 


IF(R  .LT.  O.ODO)  RETURN 
R  =  SQRT(R) 

10  IF  (L  .NE.  O.ODO)  GO  TO  15 
XJl  =  2*XJ1 
GO  TO  20 
15  V2  =  4*HZ*P 

IF  (V2  .GT.  DEP)  GO  TO  20 

XJl  =  XJ1*((0.5D0  -t-  EXP(-V2))  -I-  0.5D0) 


20  IF  (J  .GT.  1)  GO  TO  25 

CALL  CIRCV(R/S1,D,J,P,IER) 
GO  TO  30 

25  CALL  PKILL(R,S1,S2,H,K,P) 
30  FN2  =  XJ1*P 
RETURN 
END 
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SUBROUTINE  HSEXP  (X  ,E,  Y) 


EVALUATION  OF  (EXP(X)  -  1)/X 
Y  =  (EXP(X)  -  1)/X 

E  IS  AN  INPUT/OUTPUT  VARIABLE.  IF  E  IS  .GE.  0  THEN  IT  IS  ASSUMED 
THAT  E  =  EXP(X).  IN  THIS  CASE  E  IS  NOT  MODIFIED.  IF  E  IS  NEGATIVE 
THEN  E  IS  SET  TO  EXP(X)  WHEN  THIS  VALUE  IS  NEEDED  IN  HSEXP. 


DOUBLE  PRECISION  X  ,P1  ,P2  ,Q1  ,Q2  ,Q3 

DOUBLE  PRECISION  Q4  ,E  ,Y 

C - 

DATA  PI/  .914041914819518D-09/,  P2/  .238082361044469D-01/, 

♦  Q1/-.499999999085958D4-00/,  Q2/ .107141568980644D+00/, 

♦  Q3/-.119041179760821D-01/,  Q4/  .59513081 1860248D-03/ 

C - 

IF  (ABS(X)  .GT.  0.15D0)  GO  TO  5 

Y  =  ((P2»X  -h  P1)*X  +  1.0D0)/((((Q4*X  -(-  Q3)*X  -(-  Q2)*X 

♦  +  Q1)*X  -I-  l.ODO) 

RETURN 

C 

5  IF  (E  .LE.  O.DO)  E  =  EXP(X) 

IF  (X  .GT.  O.ODO)  GO  TO  10 
Y  =  ((E  -  0.5D0)  -  0.5D0)/X 
RETURN 

10  Y  =  E*(0.5D0  +  (0.5D0  -  1.0DO/E))/X 
RETURN 
END 
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SUBROUTINE  RQUAD(A,B,X11,M,N,X,Y,R3,H,K,L,S1,S2,S3,DEP,J,D) 

C - 

C  RQUAD  IS  A  QUADRATURE  ROUTINE  WHICH  USES  GAUSSIAN  MULTIPLIERS 
C  TO  OBTAIN  THE  INTEGRAL  C  FROM  A  TO  B  OF  FN2(T).  A,  B  ARE  THE 
C  LOWER  AND  UPPER  LIMITS  OF  INTEGRATION.  OUTPUT  IS  XII  =  C. 

C  INPUT  IS  M,  NO.  OF  EQUAL  SUBDIVISIONS  OF  [A,B]  WITH  SAME  ORDER 
C  GAUSSIAN  INTEGRATION  APPLIED  ON  EACH  SUBDIVISION.  M  SET  TO  1; 

C  HIGHER  VALUES  OF  M  USED  ONLY  FOR  CHECKING  PURPOSES. 

C  N,  (2N  =  ORDER  OF  GAUSSIAN  MULTIPLIERS  USED).  SET  AT  N  =  12. 

C  X(*),Y(*)-STORED  VALUES  OF  GAUSSIAN  ABSCISSAS  AND  MULTIPLIERS 
C  ON  (-1,1}.  REMAINING  ARGUMENTS  OF  THE  CALL  LINE  ARE  USED  AS 
C  INPUT  PARAMETERS  FOR  FN2,  (DEP  =  -DEPSLN(O)). 

C  J  =  0,1,2  FOR  CIRCV  OR  PKILL.  D  NEEDED  FOR  CIRCV. 


DIMENSION 

EXTERNAL 

X(l) 

FN2 

.Y{1) 

DOUBLE  PRECISION 

FN2 

,K 

,L 

DOUBLE  PRECISION 

A 

,C7 

,D 

,DEP 

,D3 

DOUBLE  PRECISION 

D4 

,E2 

,F 

,G 

,HZ 

DOUBLE  PRECISION 

R3 

,S1 

,S2 

,S3 

,T 

,TM 

DOUBLE  PRECISION 

TMl 

,TP 

,TP1 

,x 

,XI1 

,Y 

DATA  C7/.7071067811865475D0/ 

C - 

HZ  =  C7*L/S3 
5  G  =  O.DO 
K1  =  0 
D3  =  B  -  A 
D4  =  D3/M 
D3  =  D4/2 
E2  =  A  D3 
10  K1  =  Kl-t-  1 
HZ  =  C7*L/S3 

START  GAUSSIAN  INTEGRATION. 


I  =  N  -t-  I 
15  1  =  1-1 

T  =  D3*X(I) 

TMl  =  -T  E2 

TM  =  FN2(TM1,I,0,R3,H,K,L,S1,S2,S3,HZ,DEP,J,D) 
TPl  =  T  +  E2 

TP  =  FN2(TP1,I,0,R3,H,K,L,S1,S2,S3,HZ,DEP,J,D) 
F  =  Y(I)*(TM  +  TP) 

G  =  G  4-  F 
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IF  (I  .GT.  1)  GOTO  15 
E2  =  E2  -I-  D4 
IF  (K1  .NE.  M)  GOTO  10 
XII  =  D3*G 
RETURN 
END 
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SUBROUTINE  SEQHZ3(R,L,S1,S3,EPS,DEP,XK0,AA,EQS,SEQ,KKK) 


SEQ  GIVES  THE  TRIVARIATE  NORMAL  PROBABILITY  (ELLCOV)  OVER 
A  SPHERE  WITH  CENTER  (0,0, L)  AND  RADIUS  R.  THE  NORMAL  DIST- 
TRIBUTION  HAS  CORRESPONDING  STANDARD  DEVIATIONS  SI  =  S2,  S3. 
XKO  =  (L)^*2.  AA  =  EXP(-((L-R)/(SQRT(2)*MAX(S1,S3))**2)). 

EQS  =  P  FOR  EQUAL  SIGMAS  (SI  =  S3).  KKK  INDENTIFIES  PATHS, 

USED  FOR  EASE  IN  FOLLOWING  PROGRAM.  EPS  =  10*DPMPAR(1). 


SI  =  S2  AND  H  =  K  =  :. 


EXTERNAL  AERF  ,BDAW1  ,DAW  ,DXDAW  ,DXPARG 

EXTERNAL  ERF  ,ERFCR  ,ERFC0  ,HSEXP  ,EQSIG 


DOUBLE  PRECISION  AERF  ,BDAW1  ,DAW  ,DXDAW 
DOUBLE  PRECISION  DXPARG  ,ERF  ,ERFCR  ,EQSIG 


DOUBLE  PRECISION  L 


DOUBLE  PRECISION 

AA 

,BI 

,C5 

,C7 

,DEP 

,E 

DOUBLE  PRECISION 

EPS 

,EQS 

,ERL 

,ET 

,E1 

,F1 

DOUBLE  PRECISION 

R 

,s 

,SEQ 

,S0 

,S1 

,S3 

DOUBLE  PRECISION 

T 

,T0 

,T1 

,T3 

,u 

,U0 

DOUBLE  PRECISION 

U1 

,U2 

,V1 

,V2 

,V3 

,V4 

DOUBLE  PRECISION 

W9 

,x 

,XK0 

,XL1 

,XL2 

,X2 

DOUBLE  PRECISION 

Y 

,Y1 

,Y2 

,z 

,ZL 

,ZL1 

DOUBLE  PRECISION 

ZM 

,ZMN 

,ZP 

,ZPN 

,ZR 

,ZR1 

DOUBLE  PRECISION 

Z1 

DATA  B1  /.56418  95835  47756D0/,  C7/.70710  67811  86548D0/ 
DATA  C5  /11.2D0/ 


E  =  MAX  (EPS/2,1D-10) 

El  =  2DUEPS 

SEQ  =  O.DO 

ZL  =  L*C7/S3 

ZR  =  R*C7/S3 

Y1  =  4tZL*ZR 

X2  =  ZR*ZR 

VI  =  R+R/(SUS1) 


ZL  =  L/(SQRT(2)*S3)  ZR  =  R/(SQRT(2)*S3)  VI  =  R*R/(S1*SI) 


IF  (Y1  .GT.  l.Dl)  GOTO  15 
IF  (X2  .GT.  2.D0)  GOTO  15 
IF  (VI  .GT.  2.D0)  GOTO  15 
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4*ZL*ZR  .LE.  10,  ZR*ZR  .LE.  2,  VI  .LE.  2. 


FI  =  EXP(-ZUZL/2) 

S  =  .66666  66666  666667D0*F1 
T  =  S 
Y  =  2<i='l 
Z  =  .5D0*Y1 
Y1  =  Y1*F1 
Z1  =  2*X2 
J  -  0 
K  =  1 

5  K  =  K  4-  1 

Y  =  (Z*Y1  -  Z1*Y)/K 
K  =  K  -t-  1 

Y1  =  (Z*Y-  Z1*Y1)/K 
T  =  (Y/K  -  V1*T)/(K  +  2) 

S  =  S  T 

IF  (ABS(T)  .GT.  E*ABS(S))  GOTO  5 
IF  (J  .NE.  0)  GOTO  10 
J  =  1 
GOTO  5 

10  SEQ  =  B1*ZR*VUS*F1 
KKK  =  1 
RETURN 


EITHER  4*ZL*ZR  .GT.  10  .OR.  ZR*ZR  .GT.  2  .OR.  VI  .GT.  2. 


15  ZM  =  C7*(L  -  R)/S3 
KKK  =  0 

IF  (ZM  .GT.  C5)  RETURN 
V2  =  S3*S3/(S1*S1) 

XLl  =  ABS(.5D0  +  (.500  -  V2)) 

IF  (S3  .LT.  SI)  AA  =  EXP(-ZM*ZM) 
IF  (XLl  .GT.  E)  GOTO  20 
SEQ  =  EQS 


XKO  =  L*L 


KKK  =  4 

IF  (R  .GT.  L)  RETURN 

SEQ  =  EQSIG(R,XK0,L,S3,E1,DEP,AA) 

RETURN 

20  XL2  =  SQRT(XLl) 

JJJ  =  -1 

IF  (SI  .LT.  S3)  JJJ  =  1 
25  Z  =  ZR*ZR»XL1 

IF  (ZL  .NE.  O.DO  .OR.  Z  .GT.  300  )  GOTO  40 
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L  .EQ.  0  AND  R*C7/S3)**2*(ABS(1-(S3/S1)**2)).LE.  3 


K  =  3 
JJ  =  -2*JJJ 
S  =  JJ*Z/3 
T  =  S 

30  K  =  K  -I-  2 
T  =  JJ*Z*T/K 
S  =  S  -t-  T 

IF  (ABS(T)  .GT.  E*ABS(S))  GOTO  30 
IF  (S3  .GT.  SI)  GOTO  35 
EQS  =  EQSIG(R,XK0,L,S3,E1,DEP,AA) 
35  SEQ  =  EQS  -  2*ZR*B1*AA*S 
KKK  =  5 
RETURN 


L  .NE.  0  OR  R*C7/S3)*+2»(ABS(l-(S3/Sl)f*2)).GT.  3 


40  Y  =  C7*(L  -  R  +  V2*R)/S3 
Y2  =  C7*(L  +  R  -  V2*R)/S3 
IF  (V2  .LE.  l.DO)  GOTO  45 
T  =  Y 
Y  =  Y2 
Y2  =  T 

45  SO  =  Y/XL2 
FI  =  Y2/XL2 

IF  (V2  .GT.  lD-2  .OR.  V2*MAX(ZL,ZR)  .GT.  .5D0)  GOTO  65 

S3/S1  .LE.  .IDO  .AND.  MAX(ZL,ZR)*(S3/S1)**2  .LE.  1/2. 

X  =  (R-  L)*(R  L)/(2*S1*S1) 

TO  =  O.DO 

IF  (ZM  .GT.  10.5D0)  GOTO  60 
TO  =  AERF(ZL,ZR) 

ET  =  EXP(-X) 

IF  (ET  .EQ.  O.DO)  GOTO  60 


HSEXP  COMPUTES  (EXP(X)-1)/X 


CALL  HSEXP(-X,ET,SEQ) 
SEQ  =  X*SEQ*T0 


AA  =  EXP(-((L-R)*C7/S3)**2) 


ERL  =  O.DO 
ZPN  =  O.DO 
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T1  =  B1*AA 
T3  =  T1 

IF  (Y1  .GT.  -DXPARG(l))  GOTO  50 
ERL  =  EXP(-Yl) 

CALL  HSEXP(-Y1,ERL,T1) 

T1  =  T3*Y1*T1 
ZPN  =  l.DO 
50  UO  =  l.DO 
V3  =  2*V2 
U1  =  -V3*ZL 
S  =  TUUl 
ZP  =  ZL  +  ZR 
ZPN  =  l.DO 
ZMN  =  l.DO 
N  =  1 

C- - 

55  N  =  N  -t-  1 

IF  (N  .GT.  20)  GOTO  65 
ZPN  =  ZPN+ZP/N 
ZMN  =  ZMN*ZM/N 

TO  =  .5D0*T0/N  -I-  T3*(ZMN  -  ERL*ZPN) 

U2  =  V3*(-ZL*U1  (N  -  1)*U0) 

V4  =  T0*U2 

UO  =  U1 

U1  =  U2 

N  =  N  1 

ZPN  =  ZPN*ZP/N 

ZMN  =  ZMN*ZM/N 

T1  =  .5D0*T1/N  -I-  T3*(ZMN  -  ERL+ZPN) 

U2  =  V3*(-ZL*U1  -I-  (N  -  1)*U0) 

V4  =  V4  4-  U2*T1 
S  =  S  -I-  V4 
UO  =  U1 
U1  =  U2 

IF  (ABS(V4)  .GT.  ABS(SEQ  -  ET*S)*E)  GOTO  55 
SEQ  =  (SEQ  -  ET*S)/2 
KKK  =  2 
RETURN 
60  SEQ  =  TO/2 
KKK  =  3 
RETURN 


V2=(S3/S1)**2  .GT.  lD-2  .OR.  V2*MAX(ZL,ZR)  .GT.  .5D0) 


65  IF  (SO  .LT.  5D0)  GOTO  85 
Z1  =  BUZR*XL1 
W9  =  EXP(-Yl) 
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CALL  HSEXP(-Y1,W9,T) 

Z  =  4*Z*T  -  JJJ*(1  +  W9) 

U  =  AA/(2*Y*Y2) 

IF  (R  .GT.  L  .AND.  S3  .GT.  SI)  GOTO  70 
EQS  =  EQSIG(R,XK0,L,S3,E1,DEP,AA) 

70  IF  (JJJ  .GT.  0)  GOTO  80 
ZLl  =  O.DO 

IF  (W9  .EQ.  O.DO)  GOTO  75 
ZLl  =  Y*ERFCR(F1) 

75  SEQ  =  EQS  -  U*(Z1*Z  -  Y2*ERFCR(S0)  -I-  W9*ZL1) 

KKK  =  6 
RETURN 

80  SEQ  =  EQS  -  B1*U*XL1*(ZR*Z  +  BDAW1(Y,Y2,XL2,W9)) 
KKK  =  7 
RETURN 

C - 

C  USE  (13)  -  (15)  OF  REPORT  87-27 


85  SEQ  =  AERF(ZL,ZR) 

ZLl  =  ZL/XL2 
ZRl  =  XL2*ZR 
J  =  -5 

IF  (JJJ  .GT.  0)  GOTO  no 
IF  (L  .EQ.  ODO)  GOTO  105 
90  IF  (SO  .LE.  O.DO)  GOTO  100 
W9  =  EXP(-Yl) 

Y  =  -l.DO 
Y2  =  O.DO 

CALL  ERFC0(1,S0,Y,Y1) 

IF  (W9  .EQ.  O.DO)  GOTO  95 
CALL  ERFC0(1,F1,Y,Y2) 

95  SEQ  =  .5D0+(SEQ  -  AA/XL2*(Y1  -  W9*Y2)) 

KKK  =  8 
RETURN 

100  X2  =  AERF(ZL1,ZR1) 

SEQ  =  .5D0*(SEQ  -  X2*EXP(-.5D0*(R*R  -  L*L/XLl)/(SnSl))/XL2) 

KKK  =  9 

RETURN 

105  SEQ  =  SEQ/2  -  EXP(-V1/2)*ERF(ZR*XL2)/XL2 
KKK  =  10 
RETURN 

C  SI  .LT.  S3— USE  (16)  OF  REPORT  87-27 
no  J  =  -6 

IF  (L  .EQ.  ODO)  GOTO  120 

Y  =  O.DO 
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W9  =  EXP(-Yl) 

IF  (W9  .EQ.  O.DO)  GOTO  115 
Y  =  DAW(SO) 

115  SEQ  =  0.5D0*SEQ  -  B1*AA/XL2*(DAW(F1)  -  W9*Y) 
KKK  =  11 
RETURN 

120  SEQ  =  .5D0*SEQ-  2*B1*AA*ZR*DXDAW(ZR*XL2) 
KKK  =  12 
RETURN 
END 
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SUBROUTINE  SQUAD(A,B,Z1,N,X,Y,R3,H,K,L,S1,S2,S3,EPS.DEP, 

*  J,D,PY,IL) 

C. - 

C  SQUAD  IS  A  QUADRATURE  ROUTINE  WHICH  USES  GAUSSIAN  MULTIPLIERS 
C  OF  ORDER  (24)  FOR  AN  ESTIMATE  OF  THE  INTEGRAL  Z1  FROM  A  TO  B 
C  OF  FN2(T).  A,  B  ARE  THE  LOWER  AND  UPPER  LIMITS  OF  INTEGRATION. 

C  THE  OUTPUT  ARE  A,  B,  Zl,  WHERE  A,B  MAY  HAVE  NEW  VALUES,  SUCH 
C  THAT  FN2  IS  NEGLIGIBLE  ON  (A,A(NEW)),  AND  (B(NEW),B),  THAT  IS 
C  WHEN  FN2(T)  .LT.  1E-8*MIN(MAX(G,  PY,  FA,  FB),  1/2). 

C  EPS  =  10*DPMPAR(1).  DEP  =  -  DEPSLN(O).  IL  IS  INFO  PARAMETER. 

C  J  =  0,1,2  SPECIFIES  CIRCV  OR  PKILL.  CIRCV  NEEDS  D  IN  FN2. 


DIMENSION 

X(l) 

,Y(1) 

,TM(24) 

EXTERNAL 

FN2 

DOUBLE  PRECISION 

DEP 

,FN2 

,K 

,L 

DOUBLE  PRECISION 

A 

,B1 

,C7 

.D 

,D1 

DOUBLE  PRECISION 

EPS 

,E3 

,FA 

,FB 

,G 

DOUBLE  PRECISION 

HZ 

,PY 

,RZ 

,R3 

,S1 

,S2 

DOUBLE  PRECISION 

S3 

,TM 

,TM1 

.X 

,Y 

,ZT 

DOUBLE  PRECISION 

Zl 

,Z2 

DATA  B1  /.5641895835477563D0/,  C7/. 707106781 1865475D0/ 

C - 

IL  =  0 
ZT  =  0 
FA  =  O.DO 
G  =  O.DO 
NT  =  2*N 
RZ  =  C7*R3/S3 
HZ  =  C7*L/S3 

IF  (ABS(HZ  -  B)  .GT.  EPS*HZ)  GOTO  5 
Z2  =  HZ  -  RZ 

IF  (ABS(Z2  -  A)  .LE.  EPS*ABS(Z2))  IL  =  1 
5  E3  =  (B-A)/2 
D1  =  (B  4-  A)/2 
J1  =  0 
I  =  N  -I-  1 
10  1  =  1-1 
J1  =  J1  4-  1 
Z2  =  E3*X(I) 

TMl  =  -Z2  +  D1 

TM(Jl)  =  Y{1)*FN2(TM1,I,IL,R3,H,K,L,S1,S2,S3,HZ,DEP,J,D) 
TMl  =  Z2  4-  D1 
M  =  I  4-  N 
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TM(M)  =  Y(I)*FN2(TM1,M,IL,R3,H,K,L,S1,S2,S3,HZ,DEP,J,D) 
G  =  G  +  (TM(Jl)  -t-  TM(M)) 

IF  (I  .GT.  1)  GOTO  10 

IF  (IL  .EQ.  0)  FA  =  FN2(A,I,0,R3,H,K,L,S1,S2,S3,HZ,DEP,J,D) 
FB  =  FN2(B,M,0,R3,H,K,L,S1,S2,S3,HZ,DEP,J,D) 

Z1  =  E3*B1*G 

IF  (Z1  -t-  PY  .LT.  lD-37)  RETURN 

ZT  =  1D-8*MIN(MAX(G,PY,FA,FB),.5D0) 

IF  (FA  .GT.  ZT)  GOTO  25 


ATTEMPT  TO  FIND  LARGER  A 


II  =  0 

15  II  =  II  +  1 

IF  (TM(II)  .GT.  ZT)  GOTO  20 
IF  (II  .LT.  NT)  GOTO  15 
Z1  =  ODO 
RETURN 

20  IF  (II  .EQ.  1)  GOTO  25 
I  =  II  -  N  -2 
IF  (I  .GE.  0)  I  =  I  -)-  1 
JJ  =  lABS(I) 

A  =  ISIGN(1,I)*E3*X(JJ)  +  D1 
25  IF  (FB  .GT.  ZT)  RETURN 


ATTEMPT  TO  FIND  SMALLER  B 


II  =  NT  1 
30  II  =  II  -  1 

IF  (TM(Il)  .GT.  ZT)  GOTO  35 
IF  (II  .GT.  1)  GOTO  30 
35  IF  (II  .EQ.  NT)  RETURN 
I  =  II  -  N  +  1 
IF  (I  .LE.  0)  I  =  I  -  1 
JJ  =  lABS(I) 

B  =  ISIGN(1,I)*E3*X(JJ)  -t-  D1 

RETURN 

END 
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SUBROUTINE  TQUA1(A,B,N,X,R3,H,K,L,S1,S2,S3,PY,T8,DEP,F) 

C- - 

C  TQUAl  IS  A  ROUTINE  THAT  FINDS  THE  SMALLEST  VALUE  OF  T,  SAY,C  ON 
C  (A,Tl,...,Tk,...,B)  FOR  WHICH  FN2(T)  .GE.  T8,  AND  LARGEST  VALUE 
C  SAY  Cl,  FOR  WHICH  STARTING  FROM  B  IN  DECREASING  T,  FN2(C1)  .GE. 

C  T8,  WHERE  T8  =  MAX(lD2*DPMPAR(2),lD-42). 

C  TI  =  ISIGN(1,I)*E3*X(ABS(I))  -I-  D3, 

C  WITH  X(I)  THE  GAUSSIAN  ABSCISSAS  OF  ORDER  2*N,  E3  =  (B  -  A)/2, 

C  AND  D3  =  (B  +  A)/2,  X’s  =  GAUSSIAN  ABSCISSAS  OF  0(24)  USED. 

C  IF  C=T(L-f  1)  AND  C1=T(K-1),  THE  OUTPUT  A=T(L)  AND  B=T(K)  AND 

C  F  =  FN2(T(L)).  DEP  =  -DEPSLN(O). 

C - 


DIMENSION 

EXTERNAL 

DOUBLE  PRECISION 

c _ 

X(l) 

FN2 

FN2 

,K 

,L 

DOUBLE  PRECISION 

A 

,B 

,C7 

,D 

,DEP 

,D3 

DOUBLE  PRECISION 

E3 

,F 

,HZ 

,PY 

,R3 

DOUBLE  PRECISION 

SI 

,S11 

,S2 

,S22 

,S3 

,TM 

DOUBLE  PRECISION 

TMl 

,T8 

,T9 

,x 

DATA  C7/.7071067811865475D0/ 

C - 

D  =  O.DO 
HZ  =  C7*L/S3 
Sll  =  SI 
S22  =  S2 
J  =  2 
TMl  =  A 
TM  =  O.DO 

IF  (SI  .NE.  S2)  GOTO  5 
J  =  1 

D  =  SQRT(H*H  +  K*K)/S1 
GOTO  15 

5  IF  (H  +  K  .NE.  O.DO)  GOTO  15 
J  =  0 

IF  (SI  .GT.  S2)  GOTO  10 
D  =  Sll 
Sll  =  S22 
S22  =  D 

10  D  =  S22/S11 
15  NT  =  2*N  -f  1 
E3  =  (B  -  A)/2 
D3  =  (B  -I-  A)/2 
T9  =  T8 

IF  (PY  .GT.  ODO)  T9  =  MIN(T8,PY) 
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DO  20  II  =  1,  NT 
I  =  II  -  13 

IF  (I  .EQ.  0)  GOTO  20 
A  =  TMl 
F  =  TM 
JJ=IABS(I) 

TMl  =  ISIGN(1,I)*E3*X(JJ)  -I-  D3 
TM  =  FN2(TM1,I,0,R3,H,K,L,S11,S22,S3,HZ,DEP,J,D) 
IF  (TM  .GT.  T9)  GOTO  25 
20  CONTINUE 
25  TMl  =  B 
DO  30  II  =  1,  NT 
I  =  II  -  13 

IF  (I  .EQ.  0)  GOTO  30 
B  =  TMl 
JJ  =  lABS(I) 

TMl  =  -ISIGN(1,I)*E3*X(JJ)  +  D3 
IF  (TMl  .LE.  A)  GOTO  35 

TM  =  FN2(TM1,I,0,R3,H,K,L,S11,S22,S3,HZ,DEP,J,D) 
IF  (TM  .GT.  T9)  GOTO  35 
30  CONTINUE 
35  IF  (TMl  .LE.  A)  B  =  A 
RETURN 
END 
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FORTRAN  LISTINGS  FOR  ELINV3  AND  SUPPORTING  ROUTINES 


Below  we  give  a  summary  of  the  various  subprograms  that  are  used  for  computing  R  .  These 
subprograms  are  listed  in  this  appendix,  except  for  those  which  are  contained  in  NSWCLIB  [13].  The 
NSWCLIB  routines  are  identified  below  with  a  superscript  *.  All  subprograms  are  given  in  double 
precision. 

REFERENCING  OF  ROUTINES  USED  TO  COMPUTE  R 


ELINV3 

uses: 

DPMPAR*  ELLCOV 

ELLRC  GAMINV*  SUB3 

ELLRC 

uses: 

AERF*  DEPSLN* 

DPMPAR*  FCNl 

FCNl 

uses: 

AERF* 

SUB3 

uses: 

ELLCOV 
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SUBROUTINE  ELINV3(P3,HX,HY,HZ,SX,SY,SZ,R,PXD,IND,NN) 

C - 

C  (X,Y,Z)  IS  A  POINT  IN  A  CARTESIAN  COORDINATE  SYSTEM.  ELINV3 
C  RETURNS  R,  THE  RADIUS  OF  THE  SPHERE  WITH  CENTER  (HX,HY,HZ)  WHICH 
C  HOLDS  P3  OF  THE  NORMAL  ELLIPSOIDAL  DISTRIBUTION  WITH  MEAN  (0,0,0) 

C  AND  STANDARD  DEVIATIONS  SX,  SY,  SZ  IN  X,Y,Z  DIRECTIONS,  RESPECT- 
C  IVELY.  ESTIMATES  OF  P3,  FOR  A  GIVEN  R,  P(R),  ARE  OBTAINED  FROM 
C  ELLCOV.  R  IS  GENERALLY  CORRECT  TO  AT  LEAST  6  SIGNIFICANT  DIGITS. 

C  LET  E2  =  10*DPMPAR(1).  THE  INPUT  P3  SHOULD  SATISFY  lD-20  .LE.  P3 
C  LE.  MIN(1  -  E2,  .99999999).  IF  P3  IS  IN  (0,1),  AND  THE  ABOVE 
C  INEQUALIITIES  ARE  NOT  SATISFIED,  OUTPUT  R  MAY  NOT  BE  CORRECT  TO  6 
C  SIGNIFCANT  DIGITS.  IND  =  -1,  IF  P3  .LT.  MAX(lD-40,  1E6*DPMPAR(2)), 

C  P3  .NE.  0.  IND  =  1,  IF  P3  .GT.  1  -  MAX(1D-12,  E2);  R  SET  TO  -IDIO,  IF 
C  ABS(IND)  =  1.  NN  =  THE  NO.  OF  ITERATIONS  USED;  NN  .LE.  30,  NN 
C  AVERAGES  6.  IND  =  2,  IF  NN  =  30.  LET  RH  AND  RL  DENOTE  THE 
C  CURRENT  UPPER  AND  LOWER  BOUNDS  FOR  R.  IND  =  3  IF  (RH  -  RL)  .LE. 

C  MAX(E2,  1D-14)*R  AND  ABS(P(R)  -  P3)  .GT.  MAX(E2,  1E-8)*P3;  THE 
C  PRECISION  OF  THE  MACHINE  IS  NOT  ADEQUATE  TO  REVERSE  THE  LAST 
C  INEQUALITY.  IF  IND  =  4,  THE  VALUE  OF  R  IS  ACCEPTABLE  BUT  ELLCOV 
C  CANNOT  BE  COMPUTED  WITH  SUFFICIENT  ACCURACY  TO  DETERMINE  R  WITH 
C  FULL  ACCURACY. 

C  REF:  INTEGRATION  OF  THE  TRIVARIATE  NORMAL  DISTRIBUTION  OVER  AN 
C  OFFSET  SPHERE  AND  AN  INVERSE  PROBLEM.  NSWC  TR.  87-27,  2/1988. 

C  REF:  SIGNIFICANT  DIGIT  COMPUTATION  OF  THE  ELLIPSOIDAL  COVERAGE 
C  FUNCTION  AND  ITS  INVERSE.  NAVSWC  TR  91-487. 


DIMENSION 

A6(17)  ,B6(17) 

EXTERNAL 

DPMPAR  ,ELLCOV 

,ELLRC  ,GAMINV  ,SUB3 

DOUBLE  PRECISION 

DPMPAR  ,ELLCOV 

,ELLRC 

DOUBLE  PRECISION 

A6 

,B6 

,C1 

,C3 

,C4 

,D3 

DOUBLE  PRECISION 

D4 

,D5 

,E1 

,E2 

,E3 

,E4 

DOUBLE  PRECISION 

E5 

,E6 

,E8 

,F0 

,F1 

,GX 

DOUBLE  PRECISION 

GY 

,GZ 

,HX 

,HY 

,HZ 

,PH 

DOUBLE  PRECISION 

PHD 

,PL 

,PLD 

,PX 

,PXD 

,PXDO 

DOUBLE  PRECISION 

P3 

,R 

,RH 

,RL 

,RP 

,R1 

DOUBLE  PRECISION 

R2 

,s 

,SX 

,SY 

,SZ 

,V4 

DOUBLE  PRECISION 

V5 

,wx 

,WY 

,WZ 

,W2 

,W4 

DOUBLE  PRECISION 

XKO 

,XK2 

,XM0 

,X0 

,X1 

,XGAMIN 

DATA  A6(l)/lD-30/A6(2)/lD-25/A6(3)/lD-20/A6(4)/lD-15/A6(5)/lD-10/ 
+  A6(6)/lD-8/A6(7)/5D-6/A6(8)/lD-4/A6(9)/lD-2/A6(10)/lD-l/ 

+  A6(11)/3D-1/A6(12)/0.6D0/A6(13)/0.9D0/A6(14)/0.999D0/ 

+  A6(15)/.999999D0/A6(16)/.99999999D0/A6(17)/1.0D0/ 


C-3 


ooo  o  onoo 

11  I  .It  I 


NAVSWC  TR  91-487 


DATA  B6(l)/I.56D-10/B6(2)/7.23D-9/B6(3)/3.36D-7/B6(4)/1.56D-5/ 
-I-  B6(5)/7.22D-4/B6(6)/3.36D-3/B6(7)/2.66D-2/B6(8)/7.23D-2/ 

+  B6(9)/3.39D-1/B6(10)/.765D0/B6(11)/1.1933D0/B6(12)/1.717D0/ 

-I-  B6(13)/2.5005D0/B6(14)/4.0335D0/B6(15)/5.538D0/B6(16)/ 

-I-  6.35D0/B6(17)/7.7D0/ 

C - 

DATA  C4/1D-20/C3/0.797884560802865D0/C1/.5773502691896258D0/ 


C3  =  SQRT(2/PI),  Cl  =  1/SQRT(3) 


R,RL,RH,NN,PX,XO,X1,F1,PXD,FO 
IND  =  0 
R  =  ODO 

IF  (P3  .EQ.  O.ODO)  RETURN 
PXDO  =  -lODO 

XKO  =  HX*HX  +  HY*HY  +  HZ*HZ 
XK2  =  SQRT(XKO) 

S  =  MAX(SX,SY,SZ) 

El  =  1D6*DPMPAR(2) 

E2  =  10*DPMPAR(1) 

E3  =  MAX(E2,5D-9) 

E4  =  MAX(E2,1D-14) 

E5  =  MAX(E1,C4*C4) 

E6  =  MAX(E2,1D-8)*P3 

W4  =  SX*SX  -I-  SY*SY  -f  SZ*SZ 

R  =  -IDIO 

IF  (P3  .LT.  E5)  GOTO  175 

IF  (P3  .GT.  l.ODO  -  MAX(1D-12,E2))  GOTO  180 


NN  =  0 
E8  =  E6 

IF  (P3  .LT.  0.999D0)  GOTO  5 
E8  =  MAX(E2,  I.OD-9) 

IF  (P3  .LT.  0.999999D0)  GOTO  5 
E8  =  MAX(E2,  l.D-10) 

IF  (P3  .LT.  0.99999999D0)  GOTO  5 
E8  =  MAX(E2,1D-11) 


FIRST  ESTIMATE  FOR  RMIN 


5  PL  =  P3+S/C3 

RL  =  MAX(3.0D0*P3*SX*SY*SZ/C3,  PL*PL*PL) 
IF  (P3  .GE.  0.5D0)  GOTO  10 
R  =  RL**(1.D0/3.D0) 

GOTO  15 
10  R  =  XK2 

IF  (RL  .GT.  XK0*XK2)  R=B.L**(1.D0/3.D0) 


C-4 


NAVSWC  TR  91-487 


15  RL  =  R 
PL  =  ODD 
PLD  =  -  P3 

C- - 

C  FIRST  ESTIMATE  FOR  RMAX 

C- - 

DO  20  J  =  1,  17 

IF  (P3  .LE.  A6(J))  GOTO  25 
20  CONTINUE 
25  R  =  XK2  +  B6(J)+S 
RH  =  R 
PH  =  IDO 
PHD  =  IDO  -  P3 

C - 

C  GRUBB’S  ESTIMATE  FOR  R 

C - 

KG  =  0 

XMO  =  l.ODO  +  XK0/W4 
WX  =  SX*SX/W4 
WY  =  SY*SY/W4 
WZ  =  SZ*SZ/W4 
GX  =  HX/SX 
GY  =  HY/SY 
GZ  =  HZ/SZ 

V4  =  2.0D0*(WX*WX*(1.0D0  +  2.0D0*GX*GX)  +  WY*WY*(1.0D0  +  2.0D0* 
-I-  GY*GY)  +  WZ*WZ*(1.0D0  -I-  2.0D0*GZ*GZ)) 

V5  =  8.0D0*(WX*WX*WX*(I.0D0  3.0D0*GX*GX)  +  WY^WY^WY 

+  *(1.0D0  -t-  3.0D0*GY*GY)  -t-  WZ*WZ*WZ*( l.ODO  -t-  3.0D0*GZ*GZ)) 

W2  =  0.5D0*V5/V4 
V5  =  V5*V5/(V4*V4*V4) 

C- - 

12  =  0 
D3  =  I.IDO 

IF  (P3  .LE.  0.8D0)  GOTO  30 
D3  =  1.25D0 

IF  (P3  .LE.  0.9D0)  GOTO  30 
D3  =  1.9D0 
30  D4  =  P3 
PX  =  P3 

C- - 

35  RP  =  R 

IF  (PX  .EQ.  l.DO)  PX  =  l.DO  -  E8 

CALL  GAMINV(4.D0/V5,  XGAMIN,  O.DO,  PX,  l.DO  -  PX,  lERR) 

R  =  W4*((XGAM1N  -  4.0D0/V5)*W2  -|-  XMO) 

IF  (R  .LE.  RL*RL  .OR.  R  .GE.  RH*RH)  GOTO  75 
R  =  SQRT(R) 

IF  (KG  .EQ.  0)  GOTO  40 

IF  (ABS(R  -  RP)  .LT.  E3*R)  GOTO  75 
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40  CALL  SUB3(P3,R,HX,HY,HZ,SX,SY,SZ,XK0,XK2,S 
+  ,PX,PXD,RL,PL,PLD,RH,PH,PHD,NN,KK) 

KG  =  1 

45  IF  (ABS(PXD)  .LE.  E8)  RETURN 
IF  (PXD  .LT.  O.DO)  GOTO  60 

C - 

C  PX  .GE.  P3 

C - 

IF  (12  .LT.  0)  GOTO  100 
IF  (NN  .GT.  4)  GOTO  75 
IF  (PX  .GT.  10.0D0*P3)  GOTO  75 
IF  (P3  .GT.  O.OIDO  .OR.  PX  .EQ.  l.ODO)  GOTO  50 
IF  (ABS(PXD)  .GT.  0.1D0*P3)  GOTO  50 
PX  =  P3  -  2.0D0*PXD 
D4  =  P3  -  2.0D0*PXD 
GOTO  55 
50  PX  =  D4++D3 
D4  =  PX 
55  12  =  1 
GOTO  35 

C - 

C  PX  .LT.  P3 

C - 

60  D5  =  2.0D0  -  D3 

IF  (12  .GT.  0)  GOTO  100 

IF  (NN  .GT.  6  .OR.  PX  .LT.  0.01DO*P3)  GOTO  75 
IF  (P3  .GT.  O.OIDO  .OR.  PX  .EQ.  O  ODO)  GOTO  65 
IF  (ABS(PXD)  .GT.  0.1D0*P3)  GOTO  65 
PX  =  P3  -  2.0D0*PXD 
D4  =  P3  -  2.0D0*PXD 
GOTO  70 
65  PX  =  D4**D5 
D4  =  PX 
70  12  =  -1 
GOTO  35 

C  ESTIMATE  FOR  R  USING  CIRCUMSCRIBED  CUBE.  Cl  =  1/SQRT(3). 

75  R2  =  ELLRC(HX,HY,HZ,SX,SY,SZ,P3,C1*RL,RH,E8,N6) 

80  IF  (R2  .LE.  RL)  GOTO  85 
RP  =  R2 
R  =  R2 

C - 

C  COMPUTES  ELLCOV(R)  AND  MAKES  PROPER  STORAGES 

CALL  SUB3(P3,R,HX,HY,HZ,SX,SY,SZ,XK0,XK2,S 
+  ,PX,PXD,RL,PL,PLD,RH,PH,PHD,NN,KK) 

IF  (ABS(PXD)  .LE.  E8)  RETURN 
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85  R1  =  R2/C1 

IF  (R1  .LT.  RH)  RH  =  R1 

IF  ((P3  .LE.  0.2D0)  .AND.  (R2  .EQ.  RP))  GOTO  95 
90  RP  =  R 

R  =  0.5D0*(RL  +  RH) 

IF  (PL  .EQ.  ODO  .OR.  PH  .GT.  1D4*PL)  R  =  (RH  +  3*RL)/4 
IF  (ABS(RH  -  RL)  .LE.  E4*RL)  GOTO  190 

C - 

CALL  SUB3(P3,R,HX,HY,HZ,SX,SY,SZ,XK0,XK2,S 
-t-  ,PX,PXD,RL,PL,PLD,RH,PH,PHD,NN,KK) 

95  IF  (ABS(PXD)  .LE.  E8)  RETURN 

IF  (ABS(RH  -  RL)  .GT.  .1D0*R)  GOTO  90 
IF  (ABS(R  -  RP)  .GT.  5.D-3*R 
-f-  .OR.  ABS(PH  -  PL)  .GT.  l.D-3*P3)  GOTO  105 
100  RP  =  R 

R  =  RL  -  (RL  -  RH)/(PL  -  PH)*PLD 
GOTO  110 
105  RP  =  R 

R  =  0.5D0*(RL  +  RH) 

C - 

110  CALL  SUB3(P3,R,HX,HY,HZ,SX,SY,SZ,XK0,XK2,S 
-I-  ,PX,PXD,RL,PL,PLD,RH.PH,PHD,NN,KK) 

IF  (ABS(PXD).GT.  E8)  GOTO  115 
R  =  R  -  (RL  -  RH)/(PL  -  PH)*PXD 
RETURN 

115  IF  (ABS(RH  -  RL)  .LE.  E4*R)  GOTO  190 
DO  120  J2  =  1,25 
R  =  0.5D0*(RL  -t-  RH) 

IF  (PL  .EQ.  ODO  .OR.  PH  .GT.  1D4*P3)  R  =  (RH  4-  3*RL)/4 
CALL  SUB3(P3,R,HX,HY,HZ,SX,SY,SZ,XK0,XK2,S 
-t-  ,PX,PXD,RL,PL,PLD,RH,PH,PHD,NN,KK) 

IF  (ABS(PXD).LT.  E8)  RETURN 
IF  (ABS(RH  -  RL)  .LE.  E4*R)  GOTO  190 
IF  (PL  .GT.  1D-3*PH)  GOTO  125 
120  IF  (NN  .GT.  30)  GOTO  185 

C- - 

125  IK  =  0 

IF  (PL  .EQ.  ODO)  CALL  SUB3(P3,RL,HX,HY,HZ,SX,SY,SZ.XK0,XK2,S 
+  ,PX,PXD,RL,PL,PLD,RH,PH,PHD,NN,KK) 

IF  (PH  .EQ.  IDO)  CALL  SUB3(P3,RH,HX,HY,HZ,SX,SY,SZ,XK0,XK2,S 
+  ,PX,PXD,RL,PL,PLD,RH,PH,PHD,NN,KK) 

IF  (P3  .GT.  .5D0)  GOTO  130 

C  P3  .LE.  .5 

(j - 


XI  =  RL 
XO  =  RH 
FO  =  PHD 
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FI  =  PLD 
GOTO  135 

C  P3  .GT.  .5 

130  XO  =  RL 
XI  =  RH 
FO  =  PLD 
FI  =  PHD 

135  IF  (ABS(RH  -  RL)  .GT.  5.D-3*R  .OR.  ABS(PH  -  PL)  .GT.  1D-3*P3) 
*  GOTO  145 

140  R  =  0.5D0*(RL  +  RH) 

GOTO  150 

C - 

C  MODIFIED  KING’S  PROCEDURE 

C - 

145  R  =  XI  -  F1*(X1  -  X0)/(F1  -  FO) 

IF  (IK  .NE.  0)  RETURN 

150  CALL  SUB3(P3,R,HX,HY,HZ,SX,SY,SZ,XK0,XK2,S 
+  ,PX,PXD,RL,PL,PLD,RH,PH,PHD,NN,KK) 

IF  (KK  .NE.  0)  GOTO  195 

IF  (ABS(PXD)  .LE.  E8)  IK  =  1 

IF  (RH  -  RL  .LT.  E4*R)  GOTO  190 

IF  (NN  .GE.  30)  GOTO  185 

IF  (ABS(PXD  -  PXDO)  .GT.  E8)  GOTO  155 

IF  (IK  .NE.  0)  GOTO  155 

IF  (ABS^PXD)  .GT.  1D2*E8)  GOTO  155 

RETURN 

155  PXDO  =  PXD 

IF  (PXD*F1  .GT.  O.DO)  GOTO  160 

D4  =  XI 

XI  =  XO 

XO  =  D4 

D4  =  FI 

FI  =  FO 

FO  =  D4 

C - 

160  D3  =  F1/(F1  +  PXD) 

FO  =  F0*D3 
XI  =  R 
FI  =  PXD 

165  R  =  XI  -  F1*(X1  -  X0)/(F1  -  FO) 

IF  (IK  .NE.  0)  RETURN 

CALL  SUB3(P3,R,HX,HY,HZ,SX,SY,SZ,XK0,XK2,S 
+  ,PX,PXD,RL,PL,PLD,RH,PH,PHD,NN,KK) 

IF  (KK  .NE.  0)  GOTO  195 
IF  (ABS(PXD)  .LE.  E8)  IK  =  1 
IF  (RH  -  RL  .LT.  E4*R)  GOTO  190 
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IF  (NN  .GE.  30)  GOTO  185 

IF  (ABS(PXD  -  PXDO)  .GT.  E8)  GOTO  170 

IF  (IK  .NE.  0)  GOTO  170 

IF  (ABS(PXD)  .GT.  1D2*E8)  GOTO  170 

RETURN 

170  PXDO  =  PXD 

IF  (PXD*F1  .GT.  O.DO)  GOTO  160 

XO  =  XI 

FO  =  FI 

FI  =  PXD 

XI  =  R 

GOTO  145 

C - 

C  EXITS 

C  P3  .LT.  MAX(1D-40,1D6*DPMPAR(2)),  IND  =  -1,  R  =  -lElO. 

C- - 

175  IND  =  -1 
RETURN 

C- - 

C  P3  .GT.  1  -  MAX(1D-12,10*DPMPAR(1)),  IND  =  1,  R  =  -lElO. 

C - 

180  IND  =  1 
RETURN 

C- - 

C  NN  .EQ.  30 

185  IND  =  2 
RETURN 

C - 

C  RH  -  RL  .LT.  E4*R 

C - 

190  IND  =  3 
RETURN 

C- - 

C  ACCURACY  LIMITED  IN  ELLCOV 

C - 

195  IND  =  4 
RETURN 

END 


C-9 


0  0  9 


NAVSWC  TR  91-487 


FUNCTION  ELLRC(H,HK,HL,S1,S2,S3,P3,RMIN,RMAX,E8,N6) 

C- - 

C  2*ELLRC  =  LENGTH  OF  A  SIDE  OF  A  CUBE  CENTERED  AT  (H,HK,HL)  THAT 
C  CONTAINS  THE  TRIVARIATE  NORMAL  PROBABILITY  CONTENT  P3  WITH  (0,0 
C  ,0)  MEAN  AND  STANDARD  DEVIATIONS  (S1,S2,S3).  IF  N6  =  40,  RESULT 
C  SUSPECT.  RMIN  AND  RMAX  ARE  INITIAL  UPPER  AND  LOWER  BOUNDS  FOR 
C  ELLRC.  EXIT  MADE  WHEN  LATEST  ITERATE  FOR  ELLRC,  F,  SATISFIES 
C  ABS(F  -  P3)  .LT.  E8*P3.  E8  SET  IN  ELINV3. 

C - 


EXTERNAL 

AERF 

.DEPSLN 

,DPMPAR  ,FCN1 

DOUBLE  PRECISION 

AERF 

,DEPSLN 

,DPMPAR  , ELLRC 

DOUBLE  PRECISION 

A(3) 

,RA(3) 

,T(3)  ,U(3) 

DOUBLE  PRECISION 

B1 

,c 

,E 

,EPS  ,E1 

,E2 

DOUBLE  PRECISION 

E8 

,F 

,H 

,HK  ,HL 

,P3 

DOUBLE  PRECISION 

RMAX  ,RMIN 

,R0 

,R1  ,SQ 

,SQPI 

DOUBLE  PRECISION 

SI 

,S2 

,S3 

,T1  ,T2 

,V1 

C- - 

SQ  =  1.414213562373095D0 
SQPI  =  .5618958354775629D0 


SQPI  =  1/SQRT(PI) 


N6  =  0 

E  =  DPMPAR(2) 

El  =  lD-20 
El  =  EUEl 
E  =  MAX(E,E1) 

E2  =  -DEPSLN(O) 

EPS  =  MAX(5D-7,10*DPMPAR(1)) 

ELLRC  =  1.0D99 

IF  (P3  .EQ.  I.ODO)  RETURN 

ELLRC  =  O.ODO 

IF  (P3  .EQ.  O.ODO)  RETURN 

U(l)  =  1.0D0/(SQ*S1) 

U(2)  =  1.0D0/(SQ*S2) 

U(3)  =  I.0D0/(SQ*S3) 

A(l)  =  ABS(H)*U(1) 

A(2)  =  ABS(HK)*U(2) 

A(3)  =  ABS(HL)*U(3) 

RI  =  RMAX 
RO  =  RMIN 

CALL  FCN1(R1,A,U,P3,E,RA,T,F,R0,R1,IDEL) 
IF  (ABS(F  -  P3)  .GT.  E8)  GOTO  5 
ELLRC  =  Rl 
RETURN 
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5  IF  (F  .LE.  P3)  RETURN 

CALL  FCN1(R0,A,U,P3,E,RA,T,F,R0,R1,1DEL) 

IF  (ABS(F  -  P3)  .GT.  E8)  GOTO  10 

ELLRC  =  RO 

RETURN 

10  IF  (F  .GE.  P3)  RETURN 
N6  =  0 

15  ELLRC  =  0.5D0*(R1  +  RO) 

CALL  FCN1(ELLRC,A,U,P3,E,RA,T,F,R0,R1,IDEL) 

N6  =  N6  -f-  1 

IF  (N6  .EQ.  40)  RETURN 
IF  (IDEL  .NE.  0)  GO  TO  15 
IF  (ABS(F  -  P3)  .LT.  0.1D0*P3)  GOTO  30 
GO  TO  15 

20  ELLRC  =  0.5D0*(R1  RO) 

25  CALL  FCN1(ELLRC,A,U,P3,E,RA,T,F,R0,R1,IDEL) 

N6  =  N6  +  1 

IF  (N6  .EQ.  40)  RETURN 
IF  (IDEL  .NE.  0)  GO  TO  15 
30  DO  40  J  =  1,3 

T1  =  4*A(J)*RA(J) 

C  =  IDO 

IF(T1  .GT.  E2)  GOTO  35 
C  =  IDO  -F  EXP(-Tl) 

35  T2  =  (A(J)  -  RA(J)) 

RA(J)  =  EXP(-T2*T2)*C*U{J) 

40  CONTINUE 

B1  =  SQPI*(RA(1)*T(2)*T(3)  4-  RA(2)*T(1)#T(3) 

*  -I-  RA(3)*T(1)*T(2)) 

IF  (B1  .LE.  O.ODO)  GO  TO  15 
VI  =  (F  -  P3)/B1 
ELLRC  =  ELLRC  -  VI 

IF  ((ELLRC  .LT.  RO)  .OR.  (ELLRC  .GT.  Rl))  GO  TO  20 

IF  (ABS(Vl)  .GE.  5.0D-5*ELLRC)  GOTO  25 

IF  (ABS(F  -  P3)  .LT.  1.0D-3*P3)  RETURN 

IF  (ABS(Vl)  .LE.  EPS*ELLRC)  RETURN 

GO  TO  25 

END 
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SUBROUTINE  FCN1(R,A,U,P3,E,RA,T,F,R0,R1,IDEL) 

C  FCNl  USED  IN  ELLRC,  IT  GIVES  THE  PROBABILITY  F  OF  A  SHOT  FALL- 

C  ING  UNDER  A  TRIVARIATE  NORMAL  DISTRIBUTION  IN  A  CUBE  CENTERED 

C  AT  (H,HK,HL)  WITH  SIDES  OF  LENGTH  2R.  THE  DISTRIBUTION  HAS  MEAN 
C  (0,0,0)  WITH  STANDARD  DEVIATIONS  (S1,S2,S3).  IF  F  .LT.  DPMPAR(2), 

C  IDEL  =  -1.  IF  F  .GT.  1,  IDEL  =  1,  OTHERWISE  IDEL  =  0. 

C- - 

EXTERNAL  AERF 

C- - 


DOUBLE  PRECISION 

AERF  ,RA(1) 

,A(1) 

.T(l) 

,U(1) 

DOUBLE  PRECISION 

E  ,F 

,P3 

,R 

,R0  ,R1 

,v 

c - 

F  =  IDO 
IDEL  =  0 
5  DO  10  J  =  1,3 
RA(J)  =  R*U(J) 

T(J)  =  AERF(A(J),RA(J))/2 
F  =  F*T(J) 

IF  (F  .LT.  E)  GOTO  20 
10  CONTINUE 
V  =  F  -  P3 

IF  (V  .GT.  ODO)  GOTO  15 
RO  =  R 
RETURN 
15  R1  =  R 

IF  (F  .GE.  IDO)  GO  TO  25 
RETURN 
20  IDEL  =  -1 
RO  =  R 
RETURN 
25  IDEL  =  1 
RETURN 
END 
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SUBROUTINE  SUB3{P3,R3,HX,HY,HZ,SX,SY,SZ,XK0,XK2,S 
-I-  ,PX,PXD,RL,PL,PLD,RH,PH,PHD,NN,KK) 

C - 

C  SUBS  CALLS  ELLCOV.  ELLCOV  COMPARED  WITH  P3.  PX  GIVES  VALUE 
C  OF  ELLCOV.  PXD  GIVES  PX  -  P3.  PL  =  PX,  PLD  =  PXD  AND  RL  =  R3 
C  IF  PXD  .LT.  0.  PH  =  PX,  PHD  =  PXD  AND  RH  =  R3  IF  PXD  .GT.  0.  NN 
C  GIVES  NUMBER  OF  ELLCOV  CALLS.  KK  .NE.  0  MEANS  THE  ACCURACY 
C  OF  ELLCOV  NEEDED  IS  OUT  OF  REACH. 

O - 

EXTERNAL  ELLCOV 

DOUBLE  PRECISION  ELLCOV 


DOUBLE  PRECISION 

HX 

,hy 

,HZ 

,PH 

,PHD 

,PL 

DOUBLE  PRECISION 

PLD 

,PX 

,PXD 

,P3 

,RH 

,RL 

DOUBLE  PRECISION 

R3 

,s 

,SX 

,SY 

,SZ 

,XK0 

DOUBLE  PRECISION 

XK2 

KK  =  0 

PX  =  ELLCOV(R3,HX,HY,HZ,SX,SY,SZ,XKO,XK2,S) 

IF  (PX  .GT.  O.ODO)  NN  =  NN  -Hi 

PXD  =  PX  -  P3 

IF  (PXD  .GT.  O.ODO)  GO  TO  5 

RL  =  R3 

PL  =  PX 

IF  (PXD  .LT.  PLD)  KK  =  1 
PLD  =  PXD 
RETURN 
5  RH  =  R3 
PH  =  PX 

IF  (PXD  .GT.  PHD)  KK  =  1 

PHD  =  PXD 

RETURN 

END 
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